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fee CRODUCT ION 


The economics associated with ship operations have 
necessitated an examination of the losses associated with 
the motion of an automatically steered ship in a seaway. 

Four major areas -where fuel losses occur during the 
operation of a ship have been identified [Ref. 1, 2]. These 
areas on existing steam/diesel tankers are shown below: 

* Power plant and auxiliaries 
¢ Propeller efficiency 

¢ Hull resistance 

* Steering and navigation 

An optimized autopilot design would provide effective 
meeering control with associated cost Savings due to 
reducing fuel consumption. 

An appropriate computer model which represents the ship 
1s necessary for studies leading to appropriate controller 
design. Chapter 2 introduces the development of two of these 
models. 

Chapter 3 addresses the formulation of a performance 
@miterion which represents the added drag due to steering of 
the ship. 

Using the equations of motion as a model of the ship and 
PECL ION Minimization subroutine we proceed to the 


controller design for regular seas (deterministic model for 


the seaway) in Chapter 4, and for irregular seas (nondeter- 
ministic model) Pena enon ne LUNCLIOnN Minimization 
Subroutine used was BOXPLX and was programmed by R. R. 


Hilleary of the Naval Postgraduate School Computer Center 
feet. 3]. Cerin new minimum of Any arbitrary func- 
Meo, linear or nonlinear, subject to explicit constraints 
fete variables or implicit constraints on functions of the 


wart ap les . 


vk 


Chapter 6 Bitiemecices another function “minimiZaemem 
Subroutine appropriate toapeanboard Use. 

An adaptive control, which updates the controller para 
eters when the environmental conditions or the ship's cotmese 
changes 1s Studi eqgmammenaore py 

Conclusions drawn from these experiments and recommenda- 


tions for future studies are addressed in Chapter 8. 
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ie COMeute MODES OF THE SHIP 

PeeOme ly td Oat: sou diye control problem is modelling 
the process. Thus, an appropriate computer model which 
represents the ship is necessary. The best representation of 
the ship's steering dynamics is a Taylor's series expansion 
of the force and moment relationships around a_ selected 
Steady state operating point. The equations obtained in this 
way are known as the equations of motion [Ref. 4], and the 
foomu lation in Seno mateer DiEGecnail moasrmadimcated 17 Appendix 
A. This computer program was developed by using known avail- 
able data for the SL-/7 containership and by implementation 
of the scheme in Figure 2.1 [Ref. 5]. 

In this scheme the function minimization subroutine is 
fed by the yaw error Y, and rudder angle 0 , computes the 
Pemrornmance criterion J and adjusts the controller free 
Pereomleters if order £o minimize J. / 

A second model for the ship-steering dynamics represen- 
tation is the Nomoto model. Figure 2.2 indicates the second 
and third order Nomoto transfer functions while Figure 2.3 
indicates the appropriate scheme used for obtaining these 
models from the equations of motion. Appendix A includes 
the computer program used for the Nomoto third order model 
determination. 

A yaw command 1s applied as input in the scheme in 
Figure 2.3 and the difference of the signals y and Wy ms 
mameo -the function minimization subroutine is ch atte ots 
to adjust the free parameters of the Nomoto plant in order 
to minimize the performance criterion J. 

Simulation runs indicate that the resulting Nomoto 
models are obtained with resulting Jeeotosen rg IZero. 


However , in ehis SeEqdy shire equations Ot Me eae 
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representation was adopted because the system is dynamic 
and use of the Nomoto model representation implies’ addi- 
Pmemal Computer use. On the other hand, frequency domain 
StudieS were carried out using the Nomoto representation 


Since this representation 1S easier handled. 
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LIL. AN ADEQUATE PERFOR ANGE RGR eae 


A. -CRITERION BASED ONSt Ue Pie hs tr ee 


The performance criterion which characterizes propu lsum 
losses due to steering may be shown to be that derived from 
excess power consumption per unit distance caused due to 
steering [Ref. 1, 6]. The added resistance due to steering 
can be related to the surge or thrust equation wWheremmeme 


total instantaneous surge relevant to steering is 


Ax=[m+( P/2)LAX), Ivr +1/2[ (0 /2) Ax), Jv? (3s 
.+1/2[(P /2) AX ge Wimlon 


where ="MaSS OB Gomi p 

= density of sea water 

= ship's length between perpendiculars 
= [2 

ship's water speed 

= sway velocity 

= yaw rate of ship 


= rudder angle 


Pee ye ey pe Ray: oS 
iW 


~ 


vues force coefficient due to yaw/sway (posi- 
tive) 

X oo = force coefficient due to rudder angle 

(negative) 


Ry = force coefficient due to sway 


Since the sway velocity of the ship 1s, small wees 
neglect the term which includes the square of the sway 
velocity in the previous equation. From this the mean surge 


relevant to steering may be written as 
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Basted P/2)LAN),] (u4 4, /2)e0s(sP -  ) GaraZ.) 
+ [(P/2)AX U7] ( 87/2) 


Where Wa = amplitude of sway velocity 
ty = amplitude of yaw rate 
g = amplitude of rudder angle 


Oe phase difference between sway and 
r 


yaw rate 


~ 


A performance criterion for added resistance due to 


steering may be formulated as 


r 


Selim (1/2T)/(-avr+Yu? O?)dt (3.3) 
T—>0Oo O 


Vie ier craic ‘y are Cons cants 


Accurate knowledge of the MoOmMltmeanw cocliPvelents Xue and 
X ox Pe LecqiiTrederetetme accuracy, Of Such a criterion. In 
addition the criterion itself suffers from the disadvantage 
that Sway velocity measurements are not available. 


Normalizing the last equation the performance criterion will 
be 


r 


a =lim (1/2T) (-\vr+ §2)dt tor } 
NOrmM Foo O 


where A=(2[m+( 0 /2)LAJX),}/[(P /2)X 5s eal 


Mt 
Table I aundicates the values of IN for the operating 


range of speed of the ship studied. 


ey 


TARE Es 
Weighting factor 


os 


Ship's speed (knots) 


16 2A 
23 10 3477S 
52 Dio OW 


B. CRITERION BASED ONSAPPROXIMATE Aveo eee ew 


Empirical criteria based on an approximation to added 


resistance may also be derived. A semiempirical criterion 
for measuring the relative performance was developed 
(Ref. 7], based on the assumption of small amplitude oscil- 


lations around the steady-state pivot point of the ship 
during yawing at the ship/steering system natural frequency. 
This may be extended and an alternative criterion for added 


Eesms tanec Wi ape 


J=lim (1/27) fog 6 ydt (3.5) 


T—-COo 


where 
z Nw = Sn ca 


distance ae center of gravity to pivot 


rls 
i 
i 


Ceng er 
= natural frequency (closed loop ship 


steering control) 


om © 


ship's perturbation yaw angle 


é 
The values of he a function of ship's speed are given 
by Table Il. | 


A closed loop system natural frequency GW of around 0.05 


rads per sec has the potential to attenuate the effects of 
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ES eT 
4 
Weighting factor A 


Sie. Ss d (knots \ 
ion S i - (knots ) 6 720 
23 Gee 5 | 
B2 1,680 


seaway disturbance in the range of encounter angles where 
added resistance due to steering 1S important [Ref. 6]. The 
weighting factor for the operating range of the ship is 
Shown in Table III. 


TAREE, 1 LE 
Weighting factor r 


ed (knots) ; 


Saip Ss Ss 
Se 96 
Ze 28 
3 


r 
a7 
al 
RZ 


& 
6 
3 
2 


4f~COOn 


Poet wom Ss. 15 Used aS a performance criterion for this 
evICy . tiene pro canna Ole Dut Lt 1S Convenient for 
onboard use since ship's. perturbation yaw angle WV) and 


rudder oon are measurable. 


C. WEIGHTING FACTOR STUDY | 


The weighting factor IN given by Table III used in equa- 
tion 3.5, plays an important role in terms of the optimal 
controller parameters determination. Some investigation 1S 
necessary in order to verify the accuracy of the results, 


eiem@ee the values of \ Ce abe ee eeameomaebecimilned pased on 


the assumption that’ the closed loop system's natural 
Frequency is around 0.05 rads per sec [Ref. 1]. Frequency 
domain techniques were used for this purpose. Using the 


Pemeta third order model representation of the ship and 


available controller parameters from Chapter 4 for sea state 


ik 


4, encounter frequency 1.5 rads per sec, encounter angle 


150° and ship's speed 23 knots we found that the closed loop 


bandwidth of the system is 0.04 rads per sec, as indicaped 
nee eu emo Le which is not close enough toa 0.05 radsmee 
sec. 

For the same sea conditions and ship speed, with @ie 


assumption that the closed loop natural trequency of ieee 
‘system is not 0.05 rads per sec but 0.04 rads per sec, a new 
value hes en was obtained and the frequency domain tech- 
niques result in a new bandwidth 0.035 rads per Sec =e 
indicated Tne i etecs ae 

Clearly, the values of IN gBiven by Table III and used in 
this study are not the best. Unfortunately, since thei 
hydrodynamic coefficients of the SL-/7 containership are mee 
known we can't develop the surge equation and thus it is 
still impossible to determine accurate values for time 


weighting factor 
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i eeeGcUbAnwoeAs = CONTROLLER DESIGN 

We have already defined a suitable and sufficiently 
accurate ship computer model and the system's performance 
criterion, as well. The remaining task 1s to determine a 
representation of the external disturbances imparted to the 
ship by the sea, before the system's performance in a seaway 
can be evaluated. A correct model of the seaway itself is 
essential to representative modeling of forces and moments 
exerted on the ship by it. 
| At this point we will use the regular sea model as sea 
representation. The properties of regular seas are well 
defined. The wave crests are assumed to be straight, infi- 
nitely long, parallel and equally spaced with constant wave 
height. The waves progress in a direction perpendicular to 
the crest line at a uniform velocity. However the sea is 
never regular. It is a random phenomenon where waves are 
eameinually changing in height, length and breadth [Ref. 8]. 


The forces exerted by the regular sea have the form 


F=(WgRjcos( Watt ) a) 


where UL Significant wave height 
Ri = Scerrlnp LoL ce 
(a= encounter frequency 


vi phase angle 
The correspondence between sea state and wave height is 
indicated in Table IV [Ref. 9]. . 
The exciting forces RF for different encounter frequen- 
cles and encounter angles were obtained from the sea state 


generator program [Ref. 10]. 


23 


Ab ona y 


Sea state vs Wave height 


Sea state Range f e height 


One 
(Beaufort scale) ( 
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Appendix B indicates the regular seastate formulation in 
the FORTRAN program used for obtaining the controller pawaae 
a weaor 

The controller used in the entire study has one pole-one 
zero and the form is indicated in Figure 4.1. Thaes 
controller seems to have the best performance in calm waters 


and in seaway [Ref. 5]. 





Figure 4.1 Controller Used tn this seuaE 


The optimized controller parameters and the cost J ere 
23 knots speed, sea states oe different encounter 
angles and various encounter frequencies are indicated in 
Tables V,VI,VII and VIII. 
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Studying the Tables V through VIII we can draw the 
moe hOwiIne Conclusions : 

SPO 2 particular encounter angle and encounter 

frequency the higher the sea state the higher the cost. 

* For the same sea state the cost becomes smaller for 
higher encounter frequencies. 

¢ For encounter frequency 0.2 rads per sec the maximum 
cost occurs at 60° encounter angle for all tested sea 
hieates.. 

¢ For 0.6 and 0.75 rads per sec encounter frequency the 
iesimuimecostmmoccurs at 120" encounter angle for all 
Begtednsed States | 

¢ For 1.5 rads per sec encounter frequency the maximum 
cost occurs at 90° encounter angle for all tested sea 
Sieatres . 

* For 0.4 rads per sec encounter frequency the maximum 
cost occurs at 90° encounter angle for sea states 4,6 
and 7 while at sea state 9 the maximum cost occurs at 
60° encounter angle. 

Appendix C provides the computer program necessary to 
achieve the system's response. Some typical responses are 
Meme arted im Figures 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9. 
ie S obvious that as the sea state goes to heavier seas the 
rudder and yaw perturbations become larger. 

oc ECC ereTMINemenawssaccurate the controller 
Preeemeters must be for a particular situation, leads to the 
Somelusion that high accuracy isn't required. Keeping two 
parameters fixed each time and vary the third we can see 
Mepeures 4.10, 4,11, 4.12, 4.13, 4.14) that the cost doesn't 
change appreciably in the vicinity of the actual value. 

Beeotires o-oo. 4.4 5476. 4G./7, 4.8, 4.9 indicate 
that the yaw and rudder excursions are less than 1°. This 
just seems strange, though it may be because of the opti- 


Mmesation of the filter. We tried to investigate that by 


yao 


using the optimal filter for sea state 9, encounter 


frequency 1.5 rads per sec, encounter angle 1207 and rune 
in sea state 4, keeping the same encounter frequency and 
angle. The yaw and rudder excursions, even if they became 
larger, remained less than 1°. The parameters of those two 


filters are close and the reason might be the flatness of 
the cost surface. Second attempt led to more interesting 
results. Using the same filter and run it in sea state 4 
encounter frequency 0.4 rads per sec and encounter angle 


060°, the system becomes unstable (Figures 4.15, 4.16). 
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V. IRREGULAR SEAS@ = CONTROE MR Dome 


me —_—e  S 


— 


The major characteristic of the Sea is itS irrepularnimee 
This irregularity can be described by statistical methods by 
assuming that a large number of regular (sinusoidal) waves 
having different wavelengths, directions, phases and ampli- 
tudes are Superimposed to form the randomly varying Sea. 

The presence of the irregular sea was obtained by 
coupling a sea state generator program to the FORTRAN 
program as is indicated in Appendix D. The sea state gener- 
ator program generates added mass and added inertia values 
as function of the encounter frequency and also calculates 
forces and moments imparted to the shiphull by the sea. The 
forces and moments are stored in a look up table which was 
coupled to the equations of motion. The irregular sea waves 
impinging on the ship contain the total energy density spec- 
trum composed of many frequencies and the ship responds to 
an average value of added mass~ and added inertia, while in 
the regular sea the added mass and added inertia were known 
for a given encounter frequency. We decided to use values 
for added mass and added inertia corresponded to encounter 
frequency 0./5 rads per sec, Since the energy density is 
maximum in the vicinity of this frequency [Ref. 5]. This 
frequency gave us values representative of an average value 
for added mass and added inertia. 

The controller used for this study was the control@jen 
deS@mubed Sl nse napiter mm Ch aetine ee The optimized 
controller parameters and the cost J for sea states 4, 6, 7% 
9 and 0°,- 30°, 60°, 90°, 120°, IS5S0° 0 isC" encounter sansa 
are indicated in Table IX. 

Studying the Table IX we can draw the following 


CONC useoms 
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eects eamstates so.) ,ees, the higher the sea state the 
Meter enlenNGost, fOr Every particular encounter angle. 
¢ Comparing costs for sea states 4 and 6 we discover some 
Suen, Iheweest fOr a2 Specific emecounter angle in sea 
state 4 is higher than the cost for the same encounter 
om@icgeemdmesiea Staktewon Logically, we expect higher cost 
for higher sea state. 
¢ The reason for this anomaly may be the method we used 
in order to obtain the added mass and added inertia 
values. The average, we consider, might not represent 
the actual average. 
Appendix E provides the computer program necessary to 
achieve the system's response. Some typical responses are 
mieeedted in Figures 5.1, 5.2, 5.3, 5.4. 
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ieee ee ON oUD ROUTING POR ONBOARD USE 


A. GENERAL 


Peels Memeromead earlier in Chapter 1 the funetion mini- 
fmmetlom stibroutine used for these studies was BOXPLX. This 
Seecroutine will find the minimum of any function, linear or 
Memrmear, subject to explicit constraints of the variables 
Samp laiecit constraints on functions of the variables. It 
will handle a maximum of 25 variables but can handle up to 
mwevarlables with user modification. 


The variables in BOXPLX are allowed to move within a 


feasible region (n-dimensional space, where n is the number 
of variables) defined by upper and lower bounds on their 
values. The choices -for upper and lower bounds for the 


parameters are based on an understanding of the function of 
each coefficient of the system. Experience indicates that 
while accurate selections of these bounds are not necessary, 
intelligent selection of these as well as the starting 
values (guesses) can considerably reduce the computer number 
of trials needed for solution convergence. lia eesmeoron aU es@n 
was drawn trying to obtain the controller parameters for 
Maperes VY, Vi, VEIT, VIIE in Chapter 4. The function minimiza- 
tion subroutine, when starting the minimization process with 
arbitrary chosen guesses, required more than 100 trials for 
convergence while by choosing guesses close to the optimal 
parameters required more than 50 and less than 100 trials. 
Considering that every trial lasted 600 seconds (10 minutes) 
and the function minimization Subroutine requires about 60 
Samples (trials) before telling us it had found the minimum 
Bees would mean 10 hours for the control to adjust itself. 
Mere obvious reasons, such operation is not acceptable for on 


Board use. 


aS 


For obvious reasons, Such operation iS not acceptable fotmeam 


board uSe. 


B. ATTACKING THEVE Rees 


We started to investigate ways to improve this. These 
efforts include: 
© Fanci aenoaec more efficient function minimizatiom 
Subroutine 
* Studying the flame cmon nomcoss sun ie 


¢ Reducing sampling time 


GC. (Ss0OLVING 2h eRe pene 


Switching to another function minimization subroutine we 


found that the new one (ZXMWD) suffered from the same disad- 


Vania es. | 
The experiments carried out, more than two huma@iegs 
indicated that the cost surface is really flat. The BORE 


after afew trials started to focus on the minimumeess 
before it converged, it needed more than 50 trials, evens 
the guesses were close to the optimal. The reason is the way 
BOXPLK itself tries to find the minimum of a function woes 
variables. It converges when the cost FE remains unchanged 
for 2*NV consecutive trials with accuracy 10°. An efformume 
modify this termination criterion in terms of the conseeme 
tive trials was successful. Table X indicates the comparison 
between modified and unmodified BOXPLX, for sea state 6, 
encounter frequency 1.5 rads per sec and encounter angle 
1Z20am As we can observe in Table X the cost in each case 
remains almost the same while the trials required for 
convergence are dependent on the guesses made and the termi- 
nation criterion established. 

The value of the cost is in general the summation of 


incremental contributions for eadeh integration step anaes 
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ieee 
Pees eeMegieneation oan BOXPLX 


BOXPLX Guesses frivals Termination Cost 
Criterion 
Unmodified Arbitrary 100 6 0.01146720 
Unmodified Optimal ee 6 0. OPES bw 0 
Modified Poin Oe eie asia 42 3 0.01146812 
Modified Optimal le 3 OE Oar 0 
Modified aN5e lo) e Lie etelpan', i . OPS 12216 
Modified Optimal il ce 0.01146720 


therefore dependent on the total time of the simulation. 
Meneeels Important in that the optimal gain coefficients 
arrived at in this manner are not optimal for steady-state 
performance but only for the time frame covered. This should 
be adequate provided the time frame selected is long rela- 
tive to the time required for the initial condition response 
momaie out. eee oe mene sOlewcrcd tas chosen time frame 
600 seconds [Ref. 1,2]. Of course this time period is large 
and we expect steady-state behaviour faster than 10 minutes. 
Pmiplbeation studies indicate that the ship, controlled by the 
controller described in Figure 4.1, reaches the steady-state 
Situation in less than 100 seconds. So, we can reduce the 
600 seconds time frame to 200 seconds, safely. This is very 
meereait Since now the modified function minimization 
Subroutine BOXPLKX, uses samples of 200 seconds long instead 
memovd Seconds, converges in less than 30 minutes which is 
reasonable for on board use. Since the value of the cost is 
dependent on the time frame taken we expect reduced cost for 
200 seconds samples long, in any caSe. 

As we discussed earlier the BOXPLX compares the consecu- 
five trials with accuracy kG) ne Since we do not need so big 
accuracy a second modification iS necessary. We decided to 
aimee the Gxisting accuracy to 10°. 

Table XI indicates the trials required for BOXPLX 


convergence, after the second and last modification, for sea 
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TARinge «2 
Second Moditication ime so ee. 


Guesses Trigis Termination Coc 
Qi cer )] OT! 

Arbisee ai y ey 3 CP OIG RS) o) Si ehe eh?) 

Optimal 8 5 OF 003 9RS2 77 

Ar Dates ny 5 2 0.004024354 

Optimal il Z 0. 003795247 

State 6, encounter frequency 1.5 radS per sec and encounter 
eevee IL (cl By comparing Tables X and Xi we can Seem 
further improvement for the function minimization subromeume 
convergence. The difference in the cost is due to the 
different time frames used. The modified function minimiza- 


tion subroutine BOXPLA is indieated dnp rea aia 
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ieee once CONTROL 


Pee NECESSITY OF ADAPTIVITY 


The plant, Beige Soc bei whl Chim iS Supposed to be 
controlled, is normally exposed-to a time varying environ- 
ment. The ship's speed, the wave encounter angle and the sea 
Merete are changing drastically during the seaway. Since 
changes encountered are not completely predictable an 
optimum preprogrammed time-varying Grejgiiesere) Wee ah Maw 
possible. 

If we assume-and it iS apparent from the previous 
discussion in this study-that no feasible fixed parameter 
controller provides acceptable response over the entire 
performance spectrum, it is necessary that some means is 
provided for adjusting controller parameters according to 
miemoed Comditions and ship s operational characteristics. 
Perieeave CoOmtrol is thus an effort to extend basic optimum 


control concepts to these studies. 


Pee CANDIDATE ADAPTIVE SCHEMES 


Since we are looking for 14 or 2% savings in fuel cost, 
we have to feed the system with precise information. Exacer 
knowledge of the sea state, the wave encounter angle and the 
Ship's ground speed is vital for this purpose. Currently the 
Navy 1S involved ina program that will provide precision 
Mee ationm data. Garcia [|Ref. 5], provides very good infor- 
Mablon on this subject. 

An adaptive control scheme is indicated in Figure 7.1. 
Once the adaptive part of the scheme is set up the sea 
State, encounter wave angle and ship's speed are fed to the 


filters box by the appropriate sensors. iach terse box 


ae 





Controller 


(S hip) 


Seastate 
Sensor 


Decision Ship's speed 
GEV IGE Sensor 


FACOURTES 
wave angle 
Sensor 


Predetermined 


filters 





BOX PLX 


Pe Umsewa ne Adaptive Control Scheme 


ao 


Mmielides predetermined optimal filter sets for different 


G@@serete Sea States, wave angles and ship s speed. Actually, 
emose Of Lables Vig VL,VIL 


Mmeener Ss pecs ea ta lter which 


meas a look up table similar to 
and VIII. 


corresponds to discrete conditions close to those fed by the 


Piet oOutepuc Of mene 


Sensors. 
rudder angle, 
initial guesses tries 
exact sea State, 
the same time the plant is 
the optimal predetermined 
When 


tine reaches the minimum, 


mleose to the actual. 


to obtain the optimal 


wave encounter angle and ship’s speed. 


ine sumettOn minimization Subroutine accepting the 


yaw error and the predetermined filter set as 


filter for the 
ee 
eomerombledsbyetene Controller with 
set of parameters for conditions 
Ee waumeta One mms Zation, Subrou- 


eS pest me wcontrolLler with a 


new set of parameters which is the optimal set for the 


Peesele Conditions. 


But, what happens if either some or all the sensors 


provide new inputs to the system? A decision device placed 


after the sensors 


decides whether or not the change is 
meee ranle. ihis device compares the current conditions 
PoemeenOose used to obtain the controller which currently 


the system. sie 
desired percentage then a new predetermined filter is passed 


to the 


governs 1s some 


the change higher than 


SoMa ceanamre ne stunner LON mMEMimization subroutine 
tries to find the optimal controller parameters for the new 
eeeeia tion. . 

From the scheme of Figure 7.1 we can eliminate the 
predetermined filters device which provides a less expenSive 
ety 


will need a much longer time to determine the optimal filter 


system. this case the function minimization subroutine 


for rapid changes in course and speed, even if we assume 


Mie rapid changes in sea state don't occur. the 


aerial lays 


mMmectlon Minimization on 


Subset inie ml ont 


1 ae 


eemeroller set 1S obtained in every subroutine iteration 


Vici eeonturnua | ly 


imme as 1S indicated inners scheme 4 


Ieee e119 new 


and 


5 


COREG ea 






Seasfare 
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Ship's speed 
52€nsom 






Test for 


Right half 
S-plane roots Encounter 
and wave angle 
Sen sen 


tr maar ummm rene ee a ge cg pS << 


Default option 


BOX PLX 


Figure 7.2 On Line Adaptive Scheme 
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Some constraints for the controller parameters are necessary 
Mporder to avoid Operation in unstable regions. ais) the 
filter supplied by the subroutine for controlling the plant 
ict pe weested for characteristic equation roots in the 
mememt half S-plane. Goeewiae Case a derault™oction must 
exist, anda special device for such purposes 1s necesSary. 
This device will permit change in the filter parameters only 


when the new set still keeps the system in a stable situ- 


ae lon . 
Whichever adaptive scheme we adopt, we must provide for 
manual operation for the system. Manual operation will be 


desired in the following cases: 

See rivineg ports 

* Leaving ports 

* Restricted waters 

eeeavoldsecollision im sopen seas 

* Computer down 
Bee tene first two situations, Since we usually expect no 
Meemyeseas, §the optimal filter for sea state 1 seems to be 
more Suitable. Also, this filter can Serve as initial condi- 
milome tor the adaptive schemes described before, when we are 
Meavying ports. DPemgmricwmreSt Olmeune Situations the last 


operating optimal filter is the most appropriate. 


ou 


A. 


VIII. CONGLUSTONS AND REGOMME Dears ii 


CoNCLUST eh. 


The principal conclusions from this study of thew 


containership as they related to steering may be stated as 


Poul Low ce 


It is evident that a control system which providesmiaa. 
ship heading and simultaneously reduces the propulsive 
losses does exist and therefore such a controller saves 
fuel. We can't conclude how much the savings are, 
Since there is no reference for comparison between’ the 
conventional autopilot and the dattepmlctewaiere in 
addition, holds the potential for reducing the propa 
Sive losses. The literature says that savings 1% or 2% 
1s possible. 

An adaptive controller, that minimizes propulsion 
losses as ship characteristics and environmental condi- 
tions change, may be designed uSing a self-optimalizing 
technique employing a suitable performance criterion. 
Studying every particular situation we conclude that 
the cost surface is flat and therefore accurate deter- 
mination of the controller parameters 1S nor requeeaaes 
So, if we decide to use the adaptive control scheme of 
Figure 7.1we don't have to store every particular 
filter in the look up table since one filter may be 
Suitable for different ship characteristics and sea 
Gen Gd eo ae 

The weighting factor » used in the performance 
criterion equation 3.5 plays an exceptional and impor- 
tant role in the optimal controller parameters determi- 


mation. However, it is obtained from studies based on 
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MamyeisSuMPEZONS dna stherefore isn t predicted accu- 
rately. As a consequence the controller found does not 
minimize added resistance unless the proper value of 
the weighting factor r has been used. 

* The method used to obtain the average for the added 
mass and added inertia for the irregular seas studies 
might not represent the actual average. 

° The function minimization subroutine BOXPLX after two 
modifications seems to be pretty suitable for on board 


use working as a main part of the adaptive scheme. 


Pee SCOMMENDATIONS FOR FUTURE STUDIES 


The following recommendations for future work to gain a 
fuller and deeper understanding of the problem are made as 
me bOWS : 

* Some studies are necessary in order to investigate why 
PaaweeOr GhenoDEatmead Controllers in Tables V, VI, VII, 
VIII are lag and part of them are lead filters. 

* As we stated in chapter 4 the yaw and rudder excursions 
iimPeeures 6.2 through 4.9 are less than 1”. ieee 
necessary to investigate the reason for that. It might 
be because of the optimization of the filter or the 
forces and moments have not been sealed properly. 

* It is mecessary to find out the appropriate average 
values for the added mass and added inertia before we 
attempt further studies in irregular seas. 

weune Stull bydrodynamic coefficients for the SL-7 are 
necessary in order to develop and include the surge 
equation in our ship model. So far we ignore the surge 
equation and we assumed constant Smpespecd ws waile in 
reality the added resistance due to steering must 


reduce the ship speed. 


oS 


¢ By developing the surge equation in our ship model we 


may be able to determine a good value for the weighting 
factor to be used Gm the periermmamece se wee oadom 

Also, with the surge equation available in our model we 
Should be able to calculate actual energy losses and 
Savings in fuel. 

As we stated in chapter 6 the time frames (samples) 
used by the BOXPLX must be long relative to the time 
required for the initial condition response to die out. 
Reid [Ref. 1,2] has chosen time frame 600 seconds Tomes 
This time frame is long and it is necessary to finaueme 
the sufficient time frame since the time required by 
the function minimization subroutine for convergence is 


directly proportional to the sample duration. 
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APPENDIX A 
NOMOTO THIRD ORDER MODEL DETERMINATION 


//NOMOTO JOB (XXXX,XXXX),'RESEARCH' , CLASS=J 
//*MAIN ORG=NPGVM1.XXXXP 
fmsEe FORTXCG,PARM.FORT= OPT(2)' , IMSL=DP,REGION=1024K 
//FORT.SYSIN DD * 
C IN ORDER TO PERFORM SIMULATION ONLY WHEN GAINS HAVE BEEN 
C OBTAINED CHANGE XS(*) TO X(*) AND DELETE XU(*),AND XL(*) 
DIMENSION XS(4),XU(4),XL(4) 
pesne = 0 1 
po 25.13 
WS (3)=15.675 
c(i Oo 
C xXS(I) IS THE STARTING GUESS 
XL(I) IS THE LOWER LIMIT FOR THE I'TH VARIABLE 
C - XU(I) IS THE UPPER LIMIT FOR THE I'TH VARIABLE 
Sa 1) = OL 
MU(1)=1.0 
pan) = 1 <0 
2a 2 0. 
KL ( 3.) =a 
ies = 100.2 
hanno 
Mt) = 1010 0 
C A DESCRIPTION OF THE FOLLOWING PARAMETERS 
Geis DISCUSSED IN BOXPLX 
en ee oe 
la oo 
NPR=0 
NAV=0 
NV=4 
IP=0 
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C 
C 
C 


C 


LD 


30 
40 


on 


THE FOLLOWING STATEMENT MUST BE CHANGED TO 
CALL PLANT (XX) 


se 


ONLY SIMULATION IS WANTED 
CALL BOXPLX(NV,NAV,NPR,NTA,R,XS,1IP,XU, XL, YMN, IER) 
WRITE (6,25) 

FORMAT(1X,' OPTIMAL GAINS’ ,/) 
De V0 a 
WRITE(6,40)1,XS(I) 

FORMAT (1X, XC) 2” ) aaa) 
WRITE (6,77) TDIFF 

FORMAT ( 1 X@iN@@c@=' FIZ) 

STOP 

END 


SUBROUTINE PLANT(XX) SIMULATES THE SHIP 
- SUBROUTINE PLANT(XX) 


COMMON TDIFF 

REAL“8 L,L2,L3,L4,L5,L6,RXR,RXI,RYR,RYI,MZR,MZ1,RX,RY 
REAL*8 X,XDOT,Y,YDOT,U,UDOT,V, VDOT, YAW,R,RDOT,TX,TY 
REAL*8 TIME,ETIME, XUDOT, XU,XUU,XVR,XVV,XDD,WA,WE,RZ 
REAL*8 YV,YR,YD,YVVR,YVRR, YVVV, YRRR, YDDD, YVDOT , TZ 
REAL*8 NV,NR,ND,NVVR,NVRR,NVVV,NRRR,NDDD,NRDOT 

REAL*8 RHO,IZ,FX,FY,MZ,XP,MASS, DELT 

REAL*8 YAWE,YAWC,D2,D 

REAL“8 K,TP1,TP2,Z,e1,e2, Noo 

DIMENSION XX(4) 


INITIAL CONDITIGNS FOR UNI Sere wre, 
SIMULATION END LIME IN SSEConnrs 


ETIME=600. 
TIME=0.0 
ICOUNT=1 


INITIALIZE DHE COSTe FUNG are 


TDIFF=0.0 


GAIN CORFEFICEENTS 10 BEVOET Aaa 


ele) 
Z=XX(2) 
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TP1=XX(3) 
TP2=XX(4) 
C X,XDOT,Y,YDOT ARE FIX COORDINATES ON EARTH 
K=0.0 
¥=0.0 
XDOT=0.0 
YDOT=0.0 
C U,UDOT,V,VDOT ARE FIX COORDINATES ON SHIP 
V=0.0 
UDOT=0.0 
VDOT=0.0 
YAW=0.0 
R=0.0 
RDOT=0.0 
C ORDERED SPEED IN FEET/SEC 
C 38.81 FT/SEC=23 KNOTS 
UC=38.81 
C AT STEADY STATE ACTUAL SPEED (U) = COMMAND SPEED (UC) 
U=UC 
C D = RUDDER ANGLE 
D=0.0 
L=880.5 _ 
ee? 
eee 
L4sL*L3 
L5=L*L4 
L6=L*LS5 
SEA DISTURBANCE 
FORCES IN X,Y DIRECTION COMPUTED IN FORCES 
C MOMENTS IN Z- 
FX=0.0 
FY=0.0 
MZ=0 .0 
RKR=-.15744D+05 
RXI=-.19950D+06 


Oe 
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RYR=0.52365D+04 


C 
C RYI=0.18699D+06 
C MZR=-,29870D+08 
C MZI=-.37751D+07 

RXR=-.50230D+04 

RXI=0.12712D+05 

RYR=0.35290D+04 

RYI=-.31909D+05 

MZR=0.38826D+07 

MZI=-.64313D+07 
C RXR=0.28540D+04 
C RXI=-.99574D+04 
e RYR=-=.6544 10204 
C RYI=0.39595D+05 
C MZR=-.13014D+08 
C MZI=0.11348D+08 
C RXR=- .75642D+04 
C RKI=0.83497D+04 
C RYR=0.23379D+05 
C RYI=-.81502D+05 
C MZR=0 .28622D+07 
C MZI=-.19388D+08 
C RER=-.37916D+04 
C RXI=0.16381D+04 
C RYR=-.76647D+05 
C RYI=- .37685D+05 
GC MZR=- .83915D+0/7 
G MZ t= 76 D+ U7 


RX=DSQRT (RXR**2+RXI**2 ) 
RY=DSQRT(RYR**2+RYI**2 ) 
RZ=DSQRT (MZR**2+MZ1**2 ) 
TX=DATAN2 (RXI,RXR) 
TY=DATAN2 (RYI,RYR) 
TZ=DATAN2 (MZI,MZR) 
C SIGNIFICANT WAVE HEIGHT; SEA STATE I-02) 2-Cmone-2 cm 
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C 


C 


C 


oo oO Geom 


PCO Onc 10" 0.717 568-20,5,9-27.0 


WA=10.0 

PCO Meo er REOGURNG@R elma 2. -3,.4,.9,.6,.75,1.0,1.5,2.5 
WE=0.4 

FeO nem NAb CeeOhPRFIGIENTS ARE INSERTED AS PARAMETERS 
RHO=1.9876 


MASS=(.0044)*(.5*RHO*L3) 
IZ=(0.00028)*(.5*RHO*LS) 
YAWE=0.0 | 
X1=0.0 
X2=0. 
X3=0. 
K4=0. 
X5=0.0 
YAW2=0.0 
200 CONTINUE 
S=DSORT(U**2 + V*2) 
INPUT YAW COMMAND 
YAWC=0.0 
IF (TIME.GE.0.0) YAWC=(1.0/57.296) 
ERROR SIGNAL TO DRIVE RUDDER (YAW ACTUAL - YAW COMMAND) 
FOR EQUATIONS OF MOTION. 
YAWE=YAW - YAWC 
D=YAWE 


S&S 2 © 


NOMOTO 3RD ORDER PLANT 


ERROR SIGNAL TO DRIVE RUDDER (YAW COMMAND - YAW ACTUAL) 
FOR NOMOTO MODEL. 
D2=YAWC-YAW2 
olny TPL 
Ke eG eX? ) 
X4=(X3-X5)/TP2 
AXIAL FORCE HYDRODYNAMIC COEFFICIENTS (SURGE) 
KUDOT=(=. 0001 )*(.5*RHO*L3) 


C7 


XUU=(-0.0003)%*(.5*RHO*L2) 
XU=(-0.0253)*(0.5*RHO*L2*S ) 
XVR=(0.0039)*(.5*RHO*L3 ) 
XVV=(-.0012)*(.5*RHO*L2 ) 
XDD=(-0.0005 )**(.5*RHO*L2*5**2 ) 

C LATERAL FORCE HYDRODYNAMIC COEFFICIENTS (SWAY) 
YV=(-0.00758)*(.5*RHO*L2*5 ) 
YR=(0.0023)*(.5*RHO*L3*S ) 

YD= (0.00145 )*(.5*RHO*L2*S$**2) 
YVVR=(0.01)*(.5**RHO*L3/S) 
YVRR=(-0.008)*(.5*RHO*L4/S) 
YVVV=(-0.03)*(.5*RHO*L2/S) 
YRRR=(0.003)*(.5*RHO“L5/S) 
YDDD=(-0.0005)*(.5*RHO*L2*S**2 } 
YVDOT=-0.30908D+07 
YV=-0.81271D+04 

C YVDOT=-.36185D+07 

E YV=-.24757D+06 

é YVDOT=- .32890D+07 

C Weer 75 0n 

C YVDOT=- .23038D+07 

C YV=-.18267D+07 

G YVDOT=- .59800D+06 

é YV=-.13260D+07 

C MOMENT ABOUT Z-AXIS HYDRODYNAMIC COEFFICIENTS (YAW) 
NV=(-0.00213)*(.5*RHO*L3*S ) 

C NR=(-0.00105)*(.5*RHO*L4¥S ) 
ND=(-0.0007)*(.5*RHO*L3*S**2) 
NVVR=(-0.015)*(.5*RHO*L4/S) 
NVRR=(-0.008)*(.5*RHO*L5S/S) 
NVVV=(0.01)*(.5*RHO*L3/S) 
NRRR=(-0.006)*(.5*RHO*L6/S) 
NDDD=(0.0001)*(.5*RHO*L3*S**2) 

C NRDOT IS THE ADDED INERTIA TERM WHICH MUST BE CHANGED 

C FOR DIFFERENT ENCOUNTER ANGLE,SPEED,ENCOUNTER FREQUENCY 
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CC). OC) Gi) sare) 


Gi Gl Gl) Gia a) 


) 


Cas 26) 


Newer -( 90.0007 7) se eto) 
SPEED=23 KNOTS,ENCOUNTER ANGLE = 60,ENCOUNTER FREQ=0.75 
NRDOT=-.26251D+12 
NR=- .53637D+09 
New Ow =e 20125 Di 
NR=- .94970D+ 10 
NRDOT=-.18671D+12 
NR=-.46860D+11 
NRDOT=-.14518D+12 
NR=-.87538D+11 
NRDOT=-.37261D+11 
NR=-.69856D+11 
SETS SEA STATE TO ZERO 
FX=Q. 
Or 
Mz 0y 
FX=WA*RX*DCOS (WE*TIME?+TX ) 
FY=WA*RY*DCOS (WE*“TIME+TY) 
MZ=WA*RZ*DCOS (WE*TIME+TZ) 
U ACTUAL SPEED 
UC COMMANDED SPEED 
XP = PROPELLER THRUST 
XP=- XUU*UC**2 
EQUATIONS OF MOTION 
UDOT=( (XVR + MASS)*V*R + XUU*U**2 + XVV*¥V**2 
1 + XDD*D*D + FX + XP )/(MASS-XUDOT) 
VDOT=(YV*V + (YR-MASS*U)*R + YD*D + YVVR*¥V**2*R 
L + YVRR*¥V*¥R¥*2 + YVVV¥VEV%*%3 
2 + YRRR*R**3 + YDDD*D**3 + FY )/(MASS-YVDOT) 
RDOT=(NV*V + NR*¥R + ND*D + NVVR*V**2*R + NVRRYV*RY*2 
L +NVVV*V*¥*3 + NRRR*¥R*¥*3 + NDDD*D**3 + MZ )/(1IZ-NRDOT) 
WHEN TO PRINTOUT 
Piomerequnm- EO. li) "GO TO 50 
CO Lew s00 


note 


C CONVERT RADIANS TO DEGREES 


50 


YAWDEG= YAW*57.296 
RDEG=R*57.296 
RDDEG=RDOT*57 .296 
DDEG=D*57.296 
YAWC=YAWC* 57.296 
MeOUNMed: 


C TEST IF WANT TO STOP 


300 


TF (TIME.GE.EDIME GO he 3400 


C INTEGRATION SIEP Size pers 


DELT=1.0 


GC INTEGRATION 


K2=X2 ee eer 
X5=X5+X4*“DELT 
YAW2=YAW2+X5*DELT 
UU 2a) es Deed 
V=V- 7 DO lb 
R=R7 Revol pata 
YAW= YAW+R*DELT 


‘C CONVERT SHIP TO FIXED COORDINAGSee ON ee 


XDOT=U*DCOS (YAW) -V*DSIN (YAW) 
YDOT=U*DSIN(YAW)+V*DCOS (YAW) 
X=X+XDOT*DELT 

Y=Y+YDOT*DELT 

TIME=TIME+DELT 
ICOUNT=ICOUNT?+1 


C COST FUNCTION 


400 


500 


TDIFF=TDIFF+ (YAW-YAW2)** 


eO TO LOC 

CONTINUE 

WRITE(6,500) TDIFF,K,Z,TP1,TP2 

FORMAT (’ ‘qitesesCOST = esr Oe e-em 
L' Z =',F15.7,' TPl 2) oF l5e7. 2k ee ee 
RETURN 

END 


ye 


APPENDIX B 
REGULAR SEASTATE FORMULATION 


MeREGUGAINS JOB (XXXK,XXXX), 'RESEARCH’ ,CLASS=C 
//*MAIN ORG=NPGVM1L.XXXXP 
// EXEC FORTXCG,PARM.FORT='OPT(2)',IMSL=DP,REGION=1024K 
//FORT.SYSIN DD * 
C IN ORDER TO PERFORM SIMULATION ONLY WHEN GAINS HAVE BEEN 
C OBTAINED CHANGE XS(**) TO X(*) AND DELETE XU(*),AND XL(*) 
DIMENSION XS(3),XU(3),XL(3) 
XS(1)=0.9650610 
een OOgn 
XS (3)=5.6194260 
C XS(I) IS THE STARTING GUESS 
C XL(I) IS THE LOWER LIMIT FOR THE I'TH VARIABLE 
C XU(1) IS THE UPPER LIMIT FOR THE I’TH VARIABLE 
SLI) 30. 3. 
XU(1)=4.0 
nee) =O 
Muee j= 15:50 
Min 3)/=1.0 
MU (39) = 2 5.0 
C A DESCRIPTION OF THE FOLLOWING PARAMETERS 
C IS DISCUSSED IN BOXPLX 
R=9./13. 
NTA=1000 
NPR=100 
NAV=0 
NV=3 
IP=0 
C THE FOLLOWING STATEMENT MUST BE CHANGED TO 
we GALL. PLANT (X) 
C IF ONLY SIMULATION IS WANTED 
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CALL BOXPLX(NV,NAV,NPR,NTA,R,XS,IP,XU, XL, YMN, LER) 
WRITE (6,25) | 


oS FORMAT(1X,' OPTIMAL GAINS',/) 
DO 30m 

30 WRITE(6,40)1,XS(1) 

40 FORMAT (1X. "X@ 302 Dae. Hla 
STOP 
END 


SUBROUTINE PLANT(XX) 

C SUBROUTINE PLANT(XX) SIMULATES THE SHIP 
COMMON TDIFF 
REAL*8 L,L2,L3,L4,L5,L6 
REAL*8 X,XDOT,Y,YDOT,U,UDOT,V, VDOT, YAW,R,RDOT 
REAL*8 TIME,ETIME,XUDOT,XUU,XVR,XVV,XDD 
REAL*8 YV,YR,YD,YVVR,YVRR, YVVV,YRRR, YDDD , YVDOT 
REAL*8 NV,NR,ND,NVVR,NVRR,NVVV,NRRR,NDDD,NRDOT 
REAL*8 RHO,IZ,FX,FY,MZ,XP,MASS,DELT,MZI,WA,WE 
REAL“8 DYAWE,YAWE,YAWC,ISE,ISR,LAMDA,D,RYR,RYI,MZR 
REAL*8 K1,T1,T2,D,X2,DX2,S$,RX,RY,RZ,TX,TY,TZ,RXR, RXI 
DIMENSION XX(3) 


CLOSE LOOP ANALYS!IS Wither reiar 


INITIAL CONDITIONS FOR INTEGRATION 
SIMUBATION ENDER IME Nes coris 
ETIME-600 70 
pb eeald SO), (0, 
ICOUNT= 1 
C INEIETEALIZE THE COst FuUyeiien 
ISE=0.0 
ISR=0.0 
TDIFY-0-0 
LAMDA=3 ize 
C- GAEN COEFEPICIENTS TO BE soe iVi7ne 
K1=XX(1) 


Cy Gy OC) Cine 3 
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ie (2 ) 
Te 3) 
e WRITE(6,1010) K1,T1,T2 
PaO OmORMamCin, Kl = Pl507,’ Tl =',FiS.7,' T2 =',F1S.7) 
C X,XDOT,Y,YDOT ARE FIX COORDINATES ON EARTH 
X=0.0 
Y=0.0 
XDOT=0.0 
YDOT=0.0 
C U,UDOT,V,VDOT ARE FIX COORDINATES ON SHIP 
sf 10) 10 
UDOT=0.0 
VDOT=0.0 
YAW=0.0 
R=0.0 
RDOT=0.0 
C ORDERED SPEED IN FEET/SEC 
fueeo8.81 Fr/SEC=23 KNOTS 
UC=38.81 
C AT STEADY STATE ACTUAL SPEED (U) = COMMAND SPEED (UC) 
U=UC 
C D = RUDDER ANGLE 
Ce 
L=880.5 
L2=L**2 
L3=L*L*L 
L4=L*L3 
L5=L*L4 
L6=L*LS 
SEA DISTURBANCE 
FORCES IN X,Y DIRECTION COMPUTED IN FORCES 
C MOMENTS IN Z 
FX=0. 
FY=0. 
MZ=0. 


ED 


GQ) Glee Camo eG me CimeC ley TC) Ci) C) ©€) ©) €} 11) OQ ©C) QO .-).O> O©>.O 


RXR= =O 2 70s 7 Dave 
REL=0 -50Ge70- 0s 
RYR=-0.20256D+04 
RYL=O7 S077 bee 
MZR=-.14310D+08 
MZ eee 
RXR=-0.99047D+04 
RXI=.15994D+06 
RYR=- .64455D+05 
RY lL =O 6 hsv 0G 
MZR = UZ Osee ines 
MZ0==549 2040207 
RXR=- O36 76D 
RXI=0.25844D+06 
RY R= = 2.27 09 Dee 
RYL = 907 pe OG 
MZR=0.11964D+09 
MZI=0.24103D+08 
RXR=- .54639D+05 
Rki=. Z2eZ25 60000 
RYR=- .28668D+06 
RYL=0 ./7evOD aes 


MZR=O 29S 2 5 65 


MZI=0.77746D+08 
RXR=0.27268D+05 
RXI=-.71601D+05 
RYR=0.14077D+05 
RYI=-.28679D+06 
MZR=- .30892D+08 
MZI=-.53246D+08 


RY=DSORT(RYR**2+RY1L**2 ) 
RZ=DSQRT (MZR**2+MZ1**2 ) 
TX=DATAN2 (RXI, RXR) 
TY=DATAN2 (RYL,RYR) 


io 


iy Ol Gla e) 


C 
C 


TZ=DATAN2 (MZI ,MZR) 
PiCvwMure Tiley v Per TeHl ona STATE 1-0.32.2-0.75,3-2.5, 
Opes mOmG = On 7-17 s5mS-20.5. 9-27 .0 
WA=10.0 
ENCOUN TURMEREOURNCY 1,.2..3,.4,.5,.6,.75,1.0,1.5,2.5 
WE=1.5 
HYDRODYNAMIC COEFFICIENTS ARE INSERTED AS PARAMETERS 
RHO=1.9876 
Mero sree) (9 5 RHO L 3 ) 
-IZ=(0.00028)*(.5*RHO*LS ) 
YAWE=0.0 
X2=0.0 
DOr 


ZOO CONT LNUE 


S=DSQRT (U**2+V*¥*2 ) 
INPUT YAW COMMAND 

YAWC=0.0 

IF (TIME.GE.0.0) YAWC=0.0 
ERROR SIGNAL TO DRIVE RUDDER(YAW ACTUAL - YAW ORDERED) 
( COMPENSATOR FILTER ) 

YAWE=YAW - YAWC 

De Oe We 

D=K1*(T1L*DX2+X2) 
AXIAL FORCE HYDRODYNAMIC COEFFICIENTS (SURGE) 
XUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED FOR 
DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER FREQUENCY 


XUDOT=(-.0001)%*(.5*RHO*L3) 
XU= (-0.0253)*(.5*RHO*L2*S ) 
XUU=(-0.0003)*(.5*RHO*L2) 
XVR= (0.0039 )*(.5*RHO*L3 ) 
XVV=(-.0012)*(.5*RHO*L2) 
XDD= (-0.0005)*(.5*RHO*L2*S**2) 

LATERAL FORCE HYDRODYNAMIC COEFFICIENTS (SWAY) 
YV=(-0.00758 )*(.5*RHO*L2*S ) 
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Gye pee recy er Gieer 0) 6) “Gy C1ioeC) 3°) 


C2 iG). Say 3) 


YR=(0.0023)%*(.5*RHO*L3*S) 
YD=(0.00145)*(.5*RHO*L2*S**2 ) 
YVVR=(0.01)*(.5*RHO*“L3/S) 
YVRR=(-0.008)*(.5*RHO*L4/S) 
YVVV=(-0.03)*(.5*RHO*L2/S) 
YRRR=(0.003)*(.5*RHO*L5/S) 
YDDD=(-0.0005)*(.5*RHO*L2*S**2) 

YUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED FOR 

DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER FREQUENCY 


YVDOT=(-0.0039)*(.5*RHO*L3 ) 
SPEED=23 KNOTS,ENCOUNTER ANGLE=60,ENCOUNTER FREQ=0.75 
YVDOT=- .30908D+07 
YV=-.81271D+04 
YVDOT=-.36185D+07 
~YV=-.24757D+06 
YVDOT=-.32890D+07 
YV=-.11775D+07 
YVDOT2 - . 23038D+07 
YV=-.18267D+07 
YVDOT=- .59800D+06 
YV=-.13260D+07 
MOMENT ABOUT Z-AXIS HYDRODYNAMIC COEFFICIENTS (YAW) 
NV=(-0.00213)*(.5*RHO*L3%S ) 
NR=(-0.00105)*(.5**RHO*L4*S ) 
ND=(-0.0007)*(.5*RHO“L3*S***2) 
NVVR= (-0.015)*(.5*RHO*L4/S) 
NVRR=(-0.008)*(.5**RHO*L5S/S) 
NVVV=(0.01)*(.5*RHO*L3/S) 
NRRR= (-0.006)*(.5*RHO*L6/S) 
NDDD= (0.0001)**(.5**RHO*L3*S**2 ) 
NRDOT IS THE ADDED INERTIA TERM WHICH MUST BE CHANGED 
FOR DIFFERENT ENCOUNTER ANGLE,SPEED,ENCOUNTER FREQUENCY 


NRDOT=(-0.00027)*(.5*RHO*LS) 
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SPEED=23 KNOTS,ENCOUNTER ANGLE=60,ENCOUNTER FREQ=0.75 
NRDOT=-.26251D+12 
NR=- .53637D+09 
NEDO e20175D 2 
eo OO 
NRDOT=-.18671D+12 
NR=-.46860D+11 
NRDOT=-.14518D+12 
NR=-.87538D+11 
NRDOT=-.37261D+11 
NR=-.69856D+11 
REGULAR WAVE SEA STATE 
FX=WA*RX*DCOS (WE*TIME+ TX ) 
FY=WA*RY*DCOS (WE*TIME+TY ) 
MZ=WA*RZ*DCOS (WE*TIME+TZ) 
U ACTUAL SPEED 
UC COMMANDED SPEED 
XP = PROPELLER THRUST 
XP=- XUU*UC**2 
EQUATIONS OF MOTION 
UDOT=( (XVR + MASS)*V*R + XUU*U**2 + XVVEVe*2 
op beers XE) (MASS=xUDOT ) 
VDOT=(YV*V + (YR-MASS*U)*R + YD*¥D + YVVR*V*¥*2*=R 
1 + YVRR*¥V¥REED + YVVV FV w*¥ 
2 + YRRR*R**3 + YDDD*D**3 + FY )/(MASS-YVDOT) 
RDOT=(NV*¥V + NR*¥R + ND*D + NVVR*V**2*R + NVRR*®V*R**2 


WHEN TO PRINTOUT 
PRMCrcOUNT. EOL lL) Ge TO 50 
commence 

CONVERT RADIANS TO DEGREES 

50 YAWDEG= YAW*57.296 
RDEG=R*57.296 
RDDEG=RDOT*57.296 
DDHE=D--o7 296 
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YAWC=YAWC* 57.296 
IG OLONAP= IL 
C TEST IF WANT TO STOP 
300 IF (TIME.GE.ETIME) GO TO 400 
C INTEGRATION STEP SIZE DELT 
DELO 
C INTEGRATION 
U=U+UDOT*“DELT 
V=V+VDOT*“DELT 
R=R+RDOT*DELT 
YAW=YAW+R*DELT 
X2=X2+DX2*DELT 
C CONVERT SHIP TO FIXED COORDINATES ON EARTH 
G XDOT=U*DCOS ( YAW) -V*DSIN( YAW) 
C YDOT=U*DSIN (YAW) +V*DCOS (YAW) 
C X=X+XDOT*DELT 
C Y=Y+YDOT*DELT 
TIME=TIME+DELT 
ICOUNT=ICOUNT+1 


ISR=ISR + D*¥*2 
GO TO 200 
C V5=TDmgh—) COst UaUemnen 
400 TDIFF=ISE+ISR 
WRITE(6,500) ISE,ISR,TDIFF,K1,T1,T2 


500 FORMAT(’ ',1X,'ISE=',FL5.7,’ ISR=',F15.7,' TOTAL=', 
LF15.7,2K, "Ki=',F15.7,2%,'Tl=]",F15.7)20e2= oe 
RETURN 
END 


The function minimization subroutine BOX hwo eae ce 
Then the following two cards must be placed. 
//GOeSY sii be 


ale 
ea 
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APPENDIX C 
SYSTEM'S RESPONSE FOR REGULAR SEAS 


//REGURESP JOB (XXXX,XXXX), ‘RESEARCH’ ,CLASS=A 
//*MAIN ORG=NPGVML.XXXXP 
// EXEC FORTXCG,PARM.FORT='OPT(2)',IMSL=DP, REGION=1024K 
//FORT.SYSIN DD * 
C IN ORDER TO PERFORM SIMULATION ONLY WHEN GAINS HAVE BEEN 
C OBTAINED CHANGE XS(*) TO X(*) AND DELETE XU(*),AND XL(*) 
COMMON J 
DIMENSION X(3) 
X¥(1)=1.5420017 
X(2)=141.2350922 
ReneS 1094s 4 
SPCALL PLANT(X) 
C IF ONLY SIMULATION IS WANTED 
CALL PLANT(X) 
WRITE (6,25) 


25 FORMAT(1X,' OPTIMAL GAINS',/) 
DOM Outils 

30 WRITE(6,40)1I,X(L) 

40 HORMPUm@nh@meX @ 2s! = erla.7 ) 
WRITE(6,50) J 

50 HORM MNGi e@m s = 15.10) 
STOP 
END 


SUBROUTINE PLANT(XX) 
C SUBROUTINE PLANT(XX) SIMULATES THE SHIP 
COMMON TDIFF 
Rn Cues 415 , 16 
REAL*8 X,XDOT,Y,YDOT,U,UDOT,V,VDOT,YAW,R,RDOT 
REAL*8 TIME,ETIME, XUDOT,XUU,XVR,XVV,XDD 
REAL*8 YV,YR,YD,YVVR,YVRR,YVVV,YRRR,YDDD,YVDOT 


o1 


REAL*8 NV,NR,ND,NVVR,NVRR,NVVV,NRRR,NDDD,NRDOT 
REAL*8 RHO,IZ,FX,FY,MZ,XP,MASS,DELT,MZI,WA,WE 

REAL*8 DYAWE,YAWE, YAWC,ISE,ISR,LAMDA,D,RYR,RYI,MZR 
REAL“8 K1,T1,T2,D,X2,D2,S,RX,RY,RZ,TX,1Y,1Z,R eR eee 
DIMENSION XX(3) 


CLOSE LOOP ANALYSIS WilherTeieR 


INITIAL CONDITIONS FOR INTEGRATION 
SIMULATION END TIME IN SECONDS 
ETIME=600. 
TIME=0.0 
ICOUNT=1 
C INITIALIZE THE COST FUNCTION 
LSE=0.0 
ISR=0.0 
TDIFF=0.0 
LAMDA=8.128 
C GAIN COEFFICIENTS TO BE OPTIMIZED 
K1=XX(1) 
T1=XX(2) 
T2=XX(3) 
C WRITE(6,1010) Kl,T1,T2 | 
C1010 FORMAT(1X,'K1 =',F15.7,' Tl =',F15.7,’ D2 0=' | Flom 
C X,XDOT,Y,YDOT ARE FIX COORDINATES ON EARTH 
X=0.0 
Y=0.0 
XDOT=0.0 
YDOT=0.0 
C U,UDOT,V,VDOT ARE FIX COORDINATES ON SHIP 
V=0.0 — 
UpOT=0 00 
VDOT=0.0 
YAW=0.0 
R=0.0 


Cyr)" Ole) =O 


SZ 


RDOT=0.0 
C ORDERED SPEED IN FEET/SEC 
ae 63882 FIT/SEC=23 KNOTS 
WiG= 3 One Z 
C AT STEADY STATE ACTUAL SPEED (U) = COMMAND SPEED (UC) 
ENC 
C D = RUDDER ANGLE 
D=0.0 
becod 15 
ie 
L3=L*L*L 
L4=L*L3 
= ee 
on 
C SEA DISTURBANCE 
C FORCES IN X,Y DIRECTION COMPUTED IN FORCES 
C MOMENTS IN Z 
On 
Be 
MZ=0. 
RXR=-.91037D+03 
RXI=0.50896D+05 
RYR=- .20256D+04 
RYI=.18077D+06 
MZR=-.14310D+08 
MZI=-.16903D+07 
RX=DSOQRT (RXR**2+RXI**2) 
RY=DSQRT(RYR**2+RYI**2) 
RZ=DSQRT (MZR**2+MZ1I**2 ) 
TX=DATAN2 (RXI, RXR) 
TY=DATAN2 (RYI, RYR) 
TZ=DATAN2 (MZI ,MZR) 
Ge SIGNIFICANT WAVE HEIGHT; SEA STATE 1-0.32,2-0.75,3-2.5, 
fem 4=5,0,5-795.6-10.0,7-17.5,8-20.5,9-27.0 
WA=27.0 
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Oy a) a) 


C 


ENCOUNTER FREQUENCY .1,.2,.3,.4.,. Sis) um 2 ceine Onn nnSnemee 
WE=0.2 

HYDRODYNAMIC COEFFICIENTS ARE INSERTED AS PARAMETERS 
RHO=1.9876 
MASS=(.0044)*(.5*RHO*L3) 
IZ=(0.00028)*(.5*RHO*LS) 
YAWE=0.0 
X2=0.0 
DX2=0.0 


200 CONTINUE 


INPUT YAW COMMAND 

YAWC=0.0 

TF (TIMESCE JO. 0)" ye —0r 0 
ERROR SIGNAL TO DRIVE RUDDER(YAW ACTUAL - YAW ORDERED) 
( COMPENSATOR FILTER ) 

YAWE=YAW - YAWC 

Dk = (VAVE 2 \y 12 

D=Ki Tl) 
AXIAL FORCE HYDRODYNAMIC COEFFICIENTS (SURGE) 
XUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED FOR 
DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER FREQUENCY 


XUDOT=(- .0001)**(.5*RHO*L3) 
XU=(-0.0253)*(.5*RHO*L2*S ) 
XUU=(-0.0003)*(.5*RHO*L2) 
XVR= (0.0039) *(.5*RHO*L3 ) 
XVV=(-.0012)*(.5*RHO*L2 ) 
XDD= (-0.0005)*(.5*RHO*L2*S$**2 ) 
LATERAL FORCE HYDRODYNAMIC COEFFICIENTS (SWAY) 
YV=(-0.00758)*(.5*RHO*L2*S) 
YR=(0.0023)*(.5*RHO*L3*S) 
YD=(0.00145 )*(.5*RHO*L2*S%*2 ) 
YVVR=(0.01)*(.5*RHO*L3/S) 
YVRR=(-0.008)*(.5*RHO*L4/S) 
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C2 Oe ere) OC) 


Cr Gy -Oir-e 


YVVV=(-0.03)*(.5*RHO*L2/$) 
YRRR=(0.003)*(.5*RHO*L5/S) 
YDDD= (-0.0005)*(.5*RHO*L2*S**2 ) 
YUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED FOR 
DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER FREQUENCY 


YVDOT=(-0.0039)*(.5*RHO*L3 ) 
SPEED=23 KNOTS,ENCOUNTER ANGLE=60,ENCOUNTER FREQ=0.75 
YVDOT=- .30908+07 
YV=-0.81271D+04 
MOMENT ABOUT Z-AXIS HYDRODYNAMIC COEFFICIENTS (YAW) 
NV=(-0.00213)*(.5*RHO*L3*S ) 
NR=(-0.00105)*(.5*RHO*L4*s ) 
ND= (-0.0007)*(.5*RHO*L3*S**2) 
NVVR=(-0.015)*(.5*RHO*L4/S } 
NVRR=(-0.008)*(.5*RHO*LS/S) 
ey Oe ON) (oh lO Slaee 
NRRR=(-0.006)*(.5*RHO*L6/S) 
NDDD= (0.0001)*(.5*RHO*L3*S**2) 
NRDOT IS THE ADDED INERTIA TERM WHICH MUST BE CHANGED 
FOR DIFFERENT ENCOUNTER ANGLE,SPEED,ENCOUNTER FREQUENCY 


NRDOT=(-0.00027 )*(.5*RHO*LS ) 

SPEED=32 KNOTS,ENCOUNTER ANGLE=60,ENCOUNTER FREQ=0.75 
NRDOT=-0.26251D+12 
NR=-0.53637D+09 

REGULAR WAVE SEA STATE 
FX=WA*RX*DCOS (WE*TIME+TX) 
FY=WA*RY*DCOS (WE*TIME+TY) 
MZ=WA*RZ*DCOS (WE*TIME+TZ) 

U ACTUAL SPEED 

UC COMMANDED SPEED 

XP = PROPELLER THRUST 
XP=-XKUU*UC**2 

EQUATIONS OF MOTION 
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OQ) 


C 


C 


UDOT=( (XVR + MASS)*V*R + KUU*U%*2 + XVVeVr?2 
1 + XDD*D*D + FX. + XP )/(MASS-XUDOT) 
VDOT=(YV*V + (YR-MASS*U)*R + YD*D + YVVR*V**2*R 
Ll + YVRR*V*¥R**2 + YVVV*¥V33 
2 + YRRR*R**3 + YDDD*D**3 + FY )/(MASS-YVDOT) 
RDOT=(NV*V + NR*R + ND*D + NVVR*¥V**2*R + NVRRYV*R #2 
L + NVVV*V**3 + NRRR*¥R**3 + NDDD*D**3 + MZ)/(IZ-NRDOT) 
WHEN TO PRINTOUT 7 
T= (1COUNT. EO, 2)=c@ anemoe 
CO TEms00 
CONVERT RADIANS TO DEGREES 
50 YAWDEG= YAW*57.296 
RDEG=R*57.296 
RDDEG=RDOT*57.296 
DDEG=D*57.296 
YAWC=YAWC* 57.296 
WRITE (6,101) TIME, YAWDEG 
101 FORMAT(1X,F12.8,1X,F12.8) 
LCOUNT=1 
TEST IF WANT TO STOP 
300 IF (TIME.GE.ETIME) GO TO 400 
INTEGRATION STEP SIZE DELT 
DEE ae 
INTEGRATION 
U=U+UDOT*DELT 
V=V+VDOT*DELT 
R=R+RDOT“DELT 
YAW=YAW+R*DELT 
X2=X2+DX2*“DELT 
CONVERT SHIP TO FIXED COORDINATES ON EARTH 
XDOT=U*DCOS (YAW) -V*DSIN (YAW) 
YDOT=U*DSIN( YAW) +V*DCOS (YAW) 
X=X+XDOT*“DELT 
Y¥=Y¥+YDOT=DELT 
TIME=TIME+DELT 
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ICOUNT=ICOUNT+1 
ISE=ISE + LAMDA*YAWE**2 
ISR=ISR + D*¥*2 
Comme 200 
C J=TDIFF= COST FUNCTION 
400 TDIFF=ISE+ISR 
WRITE(6,500) ISE,ISR,IDIFF,K1,T1,T2 


U@MPORMAT(’ ',1X.°iSE=',F15.7,° ISR=',F15.7,° TOTAL=', 
ia nn lesen, 2X, "TL=" F15.7,2%° T2=,F15-7) 
RETURN 
END 

//GO.SYSIN DD * 


ate 
es 
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APPENDIX D 
IRREGULAR SEASTATE FORMULATION 


/ {IRREGAINS JOB (XXX, XXX) | RESEARCH /Cmugeee 

//*MAIN ORG=NPGVM1.XXXXP 

// EXEC FORTXCG, PARM.FORT='OPT(2)',IMSL=DP, REGION=1024K 
//FORT.SYSIN DD * 

IN ORDER TO PERFORM SIMULATION ONLY WHEN GAINS HAVE BEEN 
OBTAINED CHANGE XS(**) TO X(*) AND DELETE XU(*),AND XL(*) 


C 
C 


See) 
XL 
XU(I) 


DIMENSION XS(3),XU(3),XL(3) 
XS(1)=0.655751 

XS (2)=90.5483 

XS (3)=36.74847 

IS THE STARTING GUESS 

IS THE LOWER LIMIT FOR THE I'TH VARIABLE 
IS THE UPPER LIMIT FOR THE I’TH VARIABLE 
XL(1)=0.01 

Ue ees 

XL(2)=20.0 

XU(2)=180.0 

XL(3)=5.0 

XU(3)=180.0 


C A DESCRIPTION OF THE FOLLOWING PARAMETERS 
IS DISCUSSED IN BOXPLX 

R=9./13. 

NTA=1000 

NPR=100 

NAV=0 

NV=3 

hee 
THE FOLLOWING STATEMENT MUST BE CHANGED TO 
CALL PLANT(X) 
IF ONLY SIMULATION IS WANTED 


C 


C 
C 
C 
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zo 


30 
40 


C 


Cy. ©); <@)> (G2 "O) 


CALL BOXPLX(NV,NAV,NPR,NTA,R,XS,IP,XU, XL, YMN, IER) 


WRITE 


icons) 


FORMAT(1X,' OPTIMAL GAINS',/) 


bog 30 


T=1,3 


WRITE(6,40)1,XS(I) 
FORMAT(1X,'X(',12,')=',F14.7) 


Sao 
END 


SUBROUTINE PLANT(XX) 
SUBROUTINE PLANT(XX) SIMULATES THE SHIP 


COMMON 
REAL*8 
Ree ees ou 
REAL*8 
Rees 
REAL*8 
REALS 
REAL 8 
REA 


TDIFF 

Tey 164 esi 
X,XDOT,Y,YDOT,U,UDOT,V, VDOT, YAW,R,RDOT 
TIME , ETIME, XUDOT ,XUU,XVR,XVV,XDD 
YV,YR,YD,YVVR,YVRR,YVVV,YRRR,YDDD,YVDOT 
NV,NR,ND,NVVR,NVRR,NVVV,NRRR,NDDD,NRDOT 
RHO ,1Z,FX,FY,MZ,XP,MASS,DELT 

DYAWE, YAWE, YAWC,ISE,ISR,LAMDA,D 

eT ho hee? -DX2ex3 DX3.X4,CH(11),S 


DIMENSION XX(3) 


eEOse LOOP ANALYSIS WITH FILTER 


INITIAL CONDITIONS FOR INTEGRATION 
vuole D LiME IN SECONDS 
BlIME-ooe 


TIME=0. 


ILCOUNT= 


0 
IL 


PNET TALIZhetHe COsfl FUNCTION 


Sh = Ole 
ISR=0.0 


TDIFF=0.0 
LAMDA=8.128 

GAIN COEFFICIENTS TO BE OPTIMIZED 
l= Kcr) 
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ilies = Kees) 
T= oe 
6 WRITE (6, 1000) ) ke een 
C X,XDOT,Y,YDOT ARE FIX COORDINATES ON EARTH 
KeO 6 
aoe 
Domne 
DO 10 Or 
C U,UDOT,V,VDOT ARE FIX COORDINATES ON SHIP 
Ton 
UbOm= ome 
VDem—-o-0 
YAW=0.0 
R=0.0 
RDOe On 
ORDERED SPEED IN FEET/SEC 
C 38.82 FT/SEC=23 KNOTS 
UC=38. 82 
C AT STEADY STATE ACTUAL SPEED (U) = COMMAND SPEED (UC) 
ae 
C D = RUDDER ANGLE ° 
D=0.0 
L=880.5 


©) 


SEA DiSRURBANCE 
FORGES IN X,Y DIRECTION COMPUTED Tite er. 
CGC MOMENiS. LN eZ 
FX=0. 
FY=0. 
MZ=0. 
C ISEA IS A SWITCH ;ISEA=0 (CALM WATER) ISHA=1 (S54 SiAiae 


20 


C 


) Cera o©) 


C 


See =a. 
Pion Whe rGORE TINGE NT SeakeE INSERTED AS PARAMETERS 

RHO=1.98/76 

MAS o= ( wO0SGe( SRO 13 ) 

ieee revlele® Zen) “Cro + RHO“ S)) 

YAWE=0.0 

X2=0.0 
- DX2=0.0 


Z00 CONTINUE 


S=DSQRT(U**2+V**2 ) 
INPUT YAW COMMAND 

YAWC=0.0 

IF (TIME.GE.0.0) YAWC=0.0 
ERROR SIGNAL TO DRIVE RUDDER(YAW ACTUAL - YAW ORDERED) 
McONTROLLERs FILTER ) 

YAWE=YAW - YAWC 

De = VAW E22 

D=K1*(T1*DX2+X2) 


e014 FORGE HYDRODYNAMIC GOEFFICIENTS (SURGE) 


XUDOT IS THE ADDED MASS Tow Oeool BE SCHANCED FOR 
Morin ete ENCOUNTER ANGLE, SPEED | ENCOUNTER FREQUENCY 


XUDOT=(-.0001)*(.5*RHO*L3) 
KU= (-0.0253)*(.5*RHO*L2*S) 
XUU= (-0.0003)*(.5*RHO*L2) 
XVR=(0.0039)*(.5**RHO*L3 ) 
XVV=(-.0012)*(.5*RHO*L2) 
XDD=(-0.0005)*(.5*RHO“L2*5**2 ) 
LATERAL FORCE HYDRODYNAMIC COEFFICIENTS (SWAY) 
YV=(-0.00758)*(.5*RHO*L2*S) 
YR=(0.0023)*(.5*RHO*L3*S) 
YD=(0.00145)*(.5*RHO*L2*S**2 ) 
YVVR=(0.01)%*(.5*RHO*L3/S) 
YVRR=(-0.008)*(.5*RHO*L4/S) 
YVVV=(-0.03)*(.5*RHO*L2/S) 


oak 


COO) 6)" ©) 


Gea OG) wo) 


30 


YRRR=(0.003)*(.5*RHO*L5/S) 
YDDD=(-0.0005)*(.5*RHO*L2*S**2) 
YUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED FOR 
DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER FREQUENCY 


YVDOT=(-0.0039)*(.5*RHO*L3) 

SPEED=23 KNOTS, ENCOUNTER FREQUENCY =0.75 
YVDOT=-2304300.0 

MOMENT ABOUT Z-AXIS HYDRODYNAMIC COEFFICIENTS (YAW) 
NV=(-0.00213)*(.5*RHO*L3*S) 
NR=(-0.00105)*(.5*RHO*L4*S) 
ND=(-0.0007)*(.5*RHO*L3*S**2 ) 
NVVR= (-0.015)*(.5*RHO*L4/S) 
NVRR=(-0.008)*(.5*RHO*L5/S) 
NVVV=(0.01)*(.5*RHO*“L3/S) 
NRRR=(-0.006)*(.5*RHO*L6/S) 
NDDD= (0.0001)**(.5*RHO*L3*S**2 ) 

NRDOT IS THE ADDED INERTIA TERM WHICH MUST BE CHANGED 

FOR DIFFERENT ENCOUNTER ANGLE,SPEED,ENCOUNTER FREQUENCY 


NRDOT=(-0.00027)*(.5*RHO*LS ) 
SPEED=23 KNOTS, ENCOUNTER FREQUENCY =0.75 
NRDOT=-1.4518E+11 
SETS SEA STATE TO ZERO 
Th e(LSEA- EO.) GOsTomes 
FX=0. 
FY=0. 
MZ=0. 
GO TO 35 
UNIT 12 HAS THE SEA STATE DATA NAMED CH 
IT MUST BE SYNCHRONIZED BY TIME 
READ (12) CH 
FX=CH(3) 
FY=CH(4) 
MZ=CH(8) 


oe 


0} 


C 


C 


35 CONTINUE 
U ACTUAL SPEED 
UC COMMANDED SPEED 
XP = PROPELLER THRUST 
XP=-XUU*“UC**2 
EQUATIONS OF MOTION 
(PCm ea REC voR + XUUSU 2 + KV Vays? 
1 + XDD*D*D + FX + XP )/(MASS-XUDOT) | 
VDOT=(YV*V + (YR-MASS*U)*R + YD*D + YVVR*V**2"R 
Te YVRRSV=R==2 = YVVVVe"3 
2 + YRRR*¥R**3 + YDDD*D**3 + FY )/(MASS-YVDOT) 
RDOT=(NV*V + NRER + ND*D + NVVR®V%2*R + NVRREVEREH2 
1 + NVVV*V%"3 + NRRR*R**3 + NDDD*D**3 + MZ)/(IZ-NRDOT). 
WHEN TO PRINTOUT 
IF (ICOUNT.EQ.11) GO TO 50 
GO TO 300 
CONVERT RADIANS TO DEGREES 
50 YAWDEG= YAW*57.296 
RDEG=R*57.296 
RDDEG=RDOT*57.296 
DDEG=D*57.296 
YAWC=YAWC*57 .296 
LCOUNT=1 
TEST IF WANT TO STOP 
300 IF (TIME.GE.ETIME) GO TO 400 
INTEGRATION STEP SIZE DELT 
DELT=1.0 
INTEGRATION 
U=U+UDOT*DELT 
V=V+VDOT*DELT 
R=R+RDOT*DELT 
YAW=YAW+R*®DELT 
X2=X2+DK2*DELT 
CONVERT SHIP TO FIXED COORDINATES ON EARTH 
XDOT=U*DCOS (YAW) -V*DSIN( YAW) 


+ 


+ 


oa) 


YDOT=U**DSIN (YAW) +V*DCOS (YAW) 
X=X+XDOT*DELT 
Y=Y+YDOT*DELT 
TIME=TIME+DELT 
ICOUNT=ICOUNT+1 
ISE=ISE + LAMDA*YAWE***2 
ISR=ISR + D*¥*2 
eo) FO 2G 

C J=TDIFF= COST FUNCTION 
400 TDIFF=ISE+ISR 
WRITE(6,500) TDIFF,K1,T1,T2 


500 .FORMAT(‘@"?, 1X, BORAL ="SRaS.7,° Kl =u eee 
Js’ Tiles’ . FS aia TD = ee 
REWIND 12 
RETURN 
END 


The function minimization subroWsine BOXPU colle as 
Then the following three cards must be placed. 
//GO{SYSiN Bpe- 

/* | 

/ GO.FITI2Z2FOO1 DD DISP=SHRVPSNEMsseo2loGe a 
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APPENDIX E 
SYSTEM'S RESPONSE FOR IRREGULAR SEAS 


//IRRERESP JOB (XXXX,XXXX),'RESEARCH' ,CLASS=B 
/ [MAIN ORG=NPGVM1.XXXXP . 

// EXEC FRTXCLGP,IMSL=DP,REGION=1024K 
//FORT.SYSIN DD * 


C 
C 


Zo 


30 
40 


C 


Pion DER TOSPERFORM SEiMULATION ONLY WHEN GAINS HAVE BEEN 
Sel ei NED 


DIMENSION XX(3) 


ee riMal GAINS FOR CONTROLLER 


XX(1)=.9192610 
XX (2)=18.5788116 
XX (3)=9.77668 


Ios vOULENE PEANT SIMULATES THE Si=/ CONTAINERSHIP 


CALL PLANT(XX) 

WRITE(6,25) 

FORMAT(1X,'OPTIMAL GAINS" ,/) 
Porson ns 

WRITE(6,40)1,XX(L) 
HORMAMCixe MX(' 226 )=.F14.7) 
STOP 

END 


SUBROUTINE PLANT(XX) SIMULATES THE SHIP 


SUBROUTINE PLANT(XX) 

COMMON TDIFF 

REAL*8 L,L2,L3,L4,L5,L6 

REAL*™8 X,XDOT,Y,YDOT,U,UDOT,V,VDOT,YAW,R,RDOT 
REAL*8 TIME, ETIME,XUDOT,XUU,XVR,XVV,XDD 

REAL*8 YV,YR,YD,YVVR,YVRR,YVVV,YRRR,YDDD, YVDOT 
REAL*8 NV,NR,ND,NVVR,NVRR,NVVV,NRRR,NDDD,NRDOT 
REAL*8 RHO,IZ,FX,FY,MZ,XP,MASS,DELT 
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REAL*8 DYAWE,YAWE,YAWC,ISE,1SR,LAMDA,D 
REAL*8 K1,T1,T2,D,X2,DX2,8,CH(11),DX3,X3,X4 
DIMENSION XX(3) 


CLOSE LOOP GAIN Ae esr ee eee 


INITIAL CONDITIONS FOR INTEGRATION 
C SIMULATION END TIME IN SECONDS 
ETIME=600. 
TIME=0.0 
ICOUNT=1 
C INITIALIZE THE COST FUNCTION 
ISE=0.0 
ISR=0.0 
TDIFF=0.0 
LAMDA=4. 2 
C GAIN COEFFICIENTS 
K1=XX(1) 
T1=XX(2) 
T2=XX(3) 
C X,XDOT,Y,YDOT ARE FIX COORDINATES ON EARTH 
X=0.0 
Y=0.0 
XDOT=0.0 
YDOT=0.0 
C U,UDOT,V,VDOT ARE FIX COORDINATES ON SHIP 
V=0.0 
UDOT=0.0 
VDOT=0.0 
YAW=0.0 
R=0.0 
RDOT=0.0 
YAW=0.0 
C ORDERED SPEED IN FEET/SEC 
C 38.8). Bt/Ssec-23eniens 
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UC=38.81 
C AT STEADY STATE ACTUAL SPEED (U) = COMMAND SPEED (UC) 
U=UC 
C D = RUDDER ANGLE 
Demo 
L=880.5 
LD Jun? 
feel iL 
a ee 
L5=L*L4 
L6=L*L5 
C SEA DISTURBANCE 
FORCES IN X,Y DIRECTION.COMPUTED IN FORCES 
MOMENTS IN Z 
FX=0. 
FY=0, 
MZ=0. 
@eeiSEA 1S A SWITCH; ISEA=0 (CALM WATER) ISEA=1 (SEA STATE) 
ISEA=1 
C HYDRODYNAMIC COEFFICIENTS ARE INSERTED AS PARAMETERS 
RHO=1.9876 
MASS=(.0044)*(.5*RHO*L3) 
IZ=(0.00028)*(.5*RHO*LS ) 
YAWE=0.0 
X2=0.0 
DPE 
eee 
DX3=0.0 
X4=0.0 
200 CONTINUE 
S=DSQRT(U**2 + V**2) 
C INPUT YAW COMMAND 
YAWC=0.0 
IF (TIME.GE.0.0) YAWC=0.0 
C ERROR SIGNAL TO DRIVE RUDDER(YAW ACTUAL - YAW ORDERED) 


ay 


©) 


Gre) Gn Os @) 


Gh 7G). “G1 


( COMPENSATOR FILTER ) 
YAWE=YAW - YAWC 
DX2=(YAWE-X2)/T2 
X4=K1*(T1L*DX2+X2) 
DX 3= (CXGeseny ek: 
D= (T3*DX3+X3) 
AXIAL FORCE HYDRODYNAMIC COEFFICIENTS (SURGE) 


XUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED FOR 

DIFFERENT ENCOUNTER ANGLE AND SPEED. 
XUDOT=(-.0001)*(.5*RHO*L3) 

KUU= (-0.0003)**(.5*RHO*L2) 

XVR= (0.0039 )*(.5*RHO“L3) 
XVV=(-.0012)*(.5*RHO*L2 ) 

XDD= (-0.0005 )**(.5*RHO“L2*S**2 ) 

LATERAL FORCE HYDRODYNAMIC COEFFICIENTS (SWAY) 
YV=(-0.00758)**(.5*RHO“L2*S ) 
YR=(0.0023)*(.5*RHO*L3*S) 

YD= (0.00145 )**(.5*RHO*L2**S%*2 ) 

YVVR= (0.01)*(.5**RHO*L3/S) 

YVRR=(-0.008)**(.5*RHO*L4/S) 

YVVV=(-0.03)*(. 5*RHO"“L2/S) 

YRRR= (0.003)*(.5*RHO*L5/S) 

YDDD= (-0.0005)**(.5*RHO*L2**S**2) 
YUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED FOR 
DIFFERENT ENCOUNTER ANGLE AND SPEED. 


YVDOT=(-0.0039)*(.5*RHO*L3 ) 
YVDOT=-3654800.00 
MOMENT ABOUT Z-AXIS HYDRODYNAMIC COEFFICIENTS (YAW) 

NV=(-0.00213)*(.5*RHO*L3*S) 
NR=(-0.00105)**(.5*RHO*L4*S) 

ND= (-0.0007)*(.5*RHO*L3*S**2) 
NVVR=(-0.015)**(.5*RHO*L4/S) 
NVRR=(-0.008)**(.5*RHO*LS/S) 
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CHa GleO)0 


a 


O) 


30 


5D 


NVVV=(0.01)*(.5*RHO*L3/S) 
NRRR=(-0.006)*(.5*RHO*L6/S) 
NDDD=(0.0001)*(. 5*RHO*L3*S**2 ) 
NRDOT IS THE ADDED INERTIA TERM WHICH MUST BE CHANGED 
FOR DIFFERENT ENCOUNTER ANGLE AND SPEED. 


NRDOT=(-0.00027)*(.5*RHO*LS ) 
NRDOT=-2.1815E+11 
SETS SEA STATE TO. ZERO 
EP (2SEA.EO0.1) CO TO 30 
FX=0. 
FY=0. 
MZ=0. 
GO TO 35 
fe iA eA ow ioe See STATE DATA NAMED CH 
IT MUST BE SYNCHRONIZED BY TIME 
READ (12) CH 
FX= CH(3) 
FY= CH(4) 
MZ — Grice 
CONTINUE 
U ACTUAL SPEED 
UC COMMANDED SPEED 
XP = PROPELLER THRUST 
XP= - KUU*UCE*2 
EQUATIONS OF MOTION 
UDOT=( (XVR + MASS)*V*R + XUU*U%e2 + xVVeV? 
1 + XDD*D*D + FX + XP )/(MASS-XUDOT) 
VDOT=(YV*V + (YR-MASS*S)*R + YD*D + YVVREV**2"R 
1 + YVRREV*REE2 + YVVVEVI3 
1 + YRRR®=RY*3 + YDDD*D**3 + FY )/(MASS-YVDOT) 
RDOT=(NV*V + NR*¥R + ND*D + NVVR*V*¥*2*R + NVRR*¥V*¥RE*2 
1 + NVVV*V**3 + NRRR*R**3 + NDDD*D**3 + MZ)/(IZ-NRDOT) 
WHEN TO PRINTOUT 
IF (ICOUNT.EQ.2 ) GO TO 50 


S) 3) 


GO TO 300 

C CONVERT RADIANS TO DEGREES 

50 YAWDEG= YAW*57.296 
RDEG=R*57.296 
RDDEG=RDOT*57. 296 
DDEG=D*57.296 
YAWC=YAWC*57.296 
WRITE (6,100) TIME, YAWDEG 
100 FORMAT(1X,F12.8,1X,F12.8) 

ICOUNT=1 

C TEST IF WANm moMmenes 

300 IF (TIME.GE.ETIME) GO TO 400 

C INTEGRATION STEP SIZE DELT 
DELT=1. 

C INTEGRATION 
U=U+UDOT*DELT 
V=V+VDOT*DELT 
R=R+RDOT*DELT 
YAW=YAW+R*DELT 
X2=X2+DX2*DELT 
X3=X3+DX3*DELT 

C CONVERT SHIP TO FIXED COORDINATES ON EARTH 
XDOT=U*DCOS (YAW)-V*DSIN (YAW) 
YDOT=U*DSIN (YAW )+V*DCOS (YAW) 
K=X+XDOT*“DELT 
Y=Y+YDOT*DELT 
TIME=TIME+DELT 
ICOUNT=ICOUNT+1 
ISE=ISE + LAMDA*YAWE**2 
Porson aD? 
GO TO 200 

C JsTDIFF= COST FUNCTION 

400 TDIFF=ISE+ISR 
WRITE(6,500) ISE,ISR,TDIFF 

500 FORMAT('1',5X,'ISE=',F15.7,’ ISR=',F15.7, 


LOO 


nO AT = F157 ) 
STOP 
END 
//GO.SYSIN DD * - 
; oc 
fyG@- FT12F001 DD DISP=SHR,DSN=MSS.S2160.A211 


IDeyL 


APPENDIX F 
MODIFIED MINIMIZATION SUBROUTINE 


SUBROUTINE BOXPLX (NV,NAV,NPR,NTZ,RZ,XS,IP,BU,BL, YMN, IER) 
_ | 
DIMENSION V(50,50), FUN(50), SUM(25), CEN(25), XS(NV), 
LBU(NV) ,BL(NV) 


Ky = 5 

EP = 1.E-4 

NTA = 2000 

LrCNT ZC tO) en eh 

R = RZ 

IF (R. LE sOeeOr, RCE Ie ee= 1 eee 
NVT = NV+NAV 


TOMA VARSQ ae LIC i sriaticm lie me tat 
NT = 0 : | 
C CURRENT TRIAL NO. 
NPT = 0 | 
e CURRENT NO. OF PERMISSIBLE TRIALS 
NTFS = 0 | 
CURRENT NO. OF TIMES F HAS BEEN ALMOST UNCHANGED 


CHECK TEASIBREITY OF Sick ei 


CO) OO Se 


DO 4& I=1,NV 
Wee = Sn) 
LE(BlG Une. 2 Con nome 
Mat = 2 
VE) Seer el) 
GO TO 2 
lL IF (BUG) -CEavn ye coOmneme 
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©) 


Tl =e1 
Wr eG) 
2 In (NPR.CT.0) WRITE (6,49) II 
cone en a 
Cae 
IF (IP.EQ.1) GO TO 4 
BL(I) = BL(I)+AMAX1(EP,EP*ABS(BL(I))) 
BU(L) = BU(1)-AMAX1(EP,EP*“ABS(BU(I))) 
mou) = Va 


NCE = 1 
NUMBER OF CONSTRAINT EVALUATIONS 
meal 
IF (KE(V(1,1)).EQ.0) GO TO 5 
IF (NPR.LE.0) GO TO 12 
WRITE (6,50) 
EO mol? 
5 NFE = 1 


NUMBER OF VERTICES (K) = 2 TIMES NO. OF VARIABLES. 
K = (2*NV)/3 


NUMBER OF DISPLACEMENTS ALLOWED. 
NEM s= 5 NV+ LO 


Veter meOnNaEcUliyE TRIALS WITH UNCHANGED FE TO 
LOS TERMINATE 

NCT = NLIM+NV 

ALPHA = 1.3 

FK = K 

PM = shea ale, 

BETA = ALPHA+1. 


ivcUREeoeeO OF RANDOM NUMBER GENERATOR IS ODD. 


HOS 


IOR = R*1.E7 
IF (MOD(IOR,2).=0.0) IOR=IOR+101 


SET UP INITIAL VERuwens 
FUN (1). = EEG Clee 
YMN = FUN(1) 
Beoe Sale 
FUNOLD = FUN(1) 


DO Is yit-7e 
Fil =F ene 
eum y 0 

i Bel == eat ea 


END CALCULATION IF FEASIBLE GENTROID CANNOT BES POU 
TF (LIMT.GE2NEIM) GOsters 


DO 8 J=1,NV 


RANDOM NUMBER GENERATOR (RANDU) 

IQR = IQR*65539 

IF (1OR.LT.0) LOR = 1OR121G/ 4386471 

RQX = IOR 

RQX = ROQX*.4656613E-9 

V(J,I) = BL(J)}+ROX*(BU(J)-BL(J)) 

TF (1P-£0.1) V(J. 1) =41N0 C7 eee 
8 CONTINUE 


DO 10 L=1,NLIM 
NCE = NCE+l 
IF (C7 Cie ip Os Oh meer nonats 


DO 9 J=1,NV 
Vi (os  (V (S0 eerie. 
IF (LP. EORD) vr = Aaa) 
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LO 


iL. 


eZ 


iE 3 


14 


©) 


60 


15 


©) 


©) 


16 


ee = 
CONTINUE 


CONTINUE 


1 (GN Ue. 0) eter sro 

WRITE (6,51) I 

CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V,1,FUN,CEN,I) 
rm = sl 

GO TO 48 


DO 14 J=1,NV 
SUM(J) = SUM(J)+V(J,I) 
CEN(J) = SUM(J)/FI 


ieee lO wassURE PEASTBLE CENTROID FOR STARTING. 


NCE = NCE+1 

IF (KE(CEN).EQ.0) GO TO 60 
SUM(J) = SUM(J) -V(J,I) 

Om nO 

NFE = NFE+1 

FUN(1L) = FE(V(1,I1)) 
CONTINUE 


Po mOreeOr ess lt lLING OF =INITIAR COMPLEX. 


i NERO GO TO 17 
CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN,0) 


CD ineWOiks Der E Xe THE ~ J DH. 


J = 1 


On cmde=7 aK 

IF (FUN(J).GE.FUN(I)) GO TO 16 
Jeol 

CONTINUE 
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Or 1: OC) -@ 


©). -O) 


O) 


O) 


BASIC LOOP. ELIMINATE EACH WORST VERTEX IN TURN. 
IT MUST BECOME NO LONGER WORST, NOT MERELY IMPROVED. 
FIND NEXT-TO-WORST VERTEX, THE 'JN'TH ONE. 
ly cine Sal 
Pet] HO 5 eee 


po ke 
CIE 2) el Te Le 
IF (FUN(JN).GE.FUN(I)) GO TO 18 
Thal 
18 CONTINUE 


LIMT = NUMBER OF MOVES DURING THIS TRIAL TOWARD THE 
GCENTROGD DUE 10 FUNC LEO Ary ee 
iti = 


COMPUTE CENTROID AND OVER REFLECT WORST aia 


DO 19 I=1,NV 

Vile VCD) 

SUM(I) = SUM(I)-VT 

CEN(I) = SUM(1)/FKM 

VT = BETA*CEN(1I)-ALPHA*VT 

IF (IP.EQ.1) VT = AINT(VT+.5) 


INSURE THE EXPLICIT CONSTRAINTS ARE OBSERVED. 
19 V(I,J) = AMAX1(AMIN1(VT,BU(I)),BL(I)) 


NG = Nie 
CHECK FOR IMPLICIT CONSTRAINI ViGEATIiGi 


20 DOS2Z>) N= NE es 
NCE. = sNiGE 
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PE CE( VC. a) EOLO) GO TO 26 


EVERY 'KV'TH TIME, OVER-REFLECT THE OFFENDING VERTEX 
THROUGH THE BEST VERTEX. 

IF (MOD(N,KV).NE.0) GO TO 22 

CALL FBV (K,FUN,M) 


Gee =e Ny | 
VT = BETA“V(1I,M)-ALPHA*V(I,J) 
IF (IP.EQ.1) VT = AINT(VT+.5) 
21 V(I,J) = AMAX1L(AMIN1(VT,BU(I)),BL(I)) 


GO TO 24 
CONSTRAINT VIOLATION: MOVE NEW POINT TOWARD CENTROID. 


PoeDO 22) = laNny 
VT = .5*(CEN(1)+V(I,J)) 
Che ahOn lon yTs= AINE yr 5) 
ee 

23 CONTINUE 


24 NT = NT+1 
Zoe CONTINUE 


Toke 1 


CANNOT GET FEASIBLE VERTEX BY MOVING TOWARD CENTROID, 
OR BY OVER-REFLECTING THRU THE BEST VERTEX. 

i (ere) CO TON42. 

WRITE (6,52) NT,J 

CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN,J) 

GO TOr472 


een Ae OUND EVALUATE THE OBJECTIVE FUNCTION. 
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GG): “OAs O)>:61 


O) 


int GO) 


“ZO NEE 22 Nie 


FUNTRY = FE(V(1,J)) 


TEST TO-SEE IF FUNCTION VALUE HAS NOT CHANGED. 
AFO = ABS(FUNTRY-FUNOLD ) 
AMX = AMAX1(ABS(EP*FUNOLD) , EP) 


ACTIVATE THE FOLLOWING TWO STATEMENTS 
FOR DIAGNOSTIC PURPOSES ONLY. 
WRITE (6,99) J,AFO,AMX,FUNTRY,FUNOLD,FUN(J),FUN(JN), 
LNTFS,N 
99 FORMAT (1X,13,6E15.7,215) 
IF (AFO.GT.AMX) GO TO 27 
NTFS = NTFS+l 
IF (NTFS.LT.NCT) GO TO 28 
IER = 0 
IF (NPR.LE.O) GO TO 42 
WRITE (6,53) K 
CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN,0) 
GO TO 42 


2/ NTFS = Q 


IS THE NEW VERTEX NO LONGER WORST? 


28 IF (FUNTRY.LT.FUN(JN)) GO TO 34 


TRIAL VERTEX IS STILL WORST; ADJUST TOWARD CENTROID. 
EVERY 'KV'TH TIME, OVER-REFLECT THE OFFENDING VERTEX 
THROUGH THE BEST VERTEX. 

LIMT = LIMT+1 

IF (MOD(LIMT,KV).NE.0) GO TO 30 

CALL FBV (K,FUN,M) 


DO 29 I=1,NV 
VT = BETA*V(1,M)-ALPHA*V(I,J) 
IF (IP.EQ.1) VT = AINT(VT+.5) 
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OC) 1a Ci: YO) 


29 V(I,J).= AMAX1(AMIN1(VT,BU(1I)),BL(I)) 
Com Onms2 


ODO mse i= 1 Ny 
Viewer a(GEN (TeV Cl. ))) 
econ 1) Vl =mAinT(VT+.5) 
AG =n Vt 

31 CONTINUE 


Bir (LaMT.LT.NLIM) GO TO 33 


CANNOT MAKE THE 'J'TH VERTEX NO LONGER WORST BY 
DISPLACING TOWARD THE CENTROID OR BY OVER-REFLECTING 
THRU THE BEST VERTEX. 
ie 2 
Pi CiPRew bes O) CO TO 42 
WRITE (6,52) NT, J 
CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN,J) 
Ome 42 
33 NT = NT+l 
eon lo e270 


SUCCESS: WE HAVE A REPLACEMENT FOR VERTEX J. 
34 FUN(J) FUNTRY 

FUNOLD = FUNTRY 

NPT = NPT+1 


STOP AT THE 100'TH PERMISSIBLE TRIAL 
IF (MOD(NPT,100).EQ.0) GO TO 48 


DO 36 I=1,NV 
SUM(I) = 0. 


DO 35 N=1,K 


OD 


S) 


Cr Ca, “Qe 7 Ce 1047 


35 


36 


B) // 
38 


a 


SUM(I1) 


GEN GIS) 
CONTINUE 


meas 0 
GO TO 39 


DO 38 I=1,NV 


SUM(I)+V(I,N) 


SUM(1)/FK 


SUM(I) = SUM(1)+V(I,J) 


LC = J 


LF. (NEReLE Ome OmmOomcG 
IF (MOD(NPT,NPR).NE.0) GO TO 40 


CALL BOUT (NT,NPT,NFE,NCE,NV{ NVI ,V (KK ORUN ees 


HAS THE MAX. NUMBER, OF —TRLALS BEEN REACHED) Liou 


CONVERGENCE 


TF NOI, "GO TO NEW e ie 


40 Dre(NEMeos Nie) CO. moncal 


NEATH TO=“wORST VERDE GNOW ee COME Se VOns i 


J = JN 
GO TO 17 
4] TER = 3 


IF (NPR.GT.O) WRITE (6,54) 


COLLECTOR POINT FOR ALL ENDINGS. 
1) CANNOT DEVELOP FEASIBLE VERTEX. 


2) CANNOT DEVELOP A NO-LONGER-WORST VERTEX. 


3) FUNCTION 
4) LIMIT ON 


VALUE UNCHANGED FOR K TRIALS: 
TREALS  REaS De 


5) CANNOT FIND FEASEBEE VERGE oto. 
42 CONTINUE 
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das 
aK 
aps 
Jas 
Bee 


re WwW O WM - 


FIND BEST VERTEX. 
CALL FBV (K,FUN,M) 
Moa ( LER SCE. S ) GO. T0144 
RESTART IF THIS SOLUTION IS SIGNIFICANTLY BETTER THAN 
THE PREVIOUS OR IF THIS IS THE FIRST TRY. 
IF (NPR.LE.O) GO TO 43 
WRITE (6,55) (M,YMN,FUN(M)) 
43 IF (FUN(M).GE.YMN) GO TO 47 
IF (ABS(FUN(M)-YMN).LE.AMAX1(EP,EP*YMN)) GO TO 47 
GIVE IT ANOTHER TRY UNLESS LIMIT ON TRIALS REACHED. 
44 YMN = FUN(M) 
FUN(1) = FUN(M) 


DO 45 I=1,NV 


CEN(L) = V(I,M) 

SUM(I) = V(I,M) 
45 V(I,1) = V(I,M) 

DO 46 I=1,NVT 


46 XS(I) = V(I,M) 


PEM GhEwrr. 3) CO TO 6 

47 IF (NPR.LE.O) GO TO 48 
CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,V(1,M),-1) 
WRITE (6,56) FUN(M) 

48 RETURN 


49 FORMAT (SOHOINDEX AND DIRECTION OF OUTLYING 
LVARIABLE AT STARTIS) 

50 FORMAT (SOHOIMPLICIT CONSTRAINT VIOLATED AT START. 
LDEAD END.) 

51 FORMAT ('OCANNOT FIND FEASIBLE',14,'TH VERTEX OR 
LCENTROID AT START.') 

52 FORMAT (1lOHOAT TRIAL 14,54H CANNOT FIND FEASIBLE 


JE AI 


1VERTEX WHICH IS NO 
2LONGER WORST,14,15X,'RESTART FROM BEST VERTEX. ' ) 
53 FORMAT (4OHOFUNCTION HAS BEEN ALMOST UNCHANGED 
1FOR I5,7H TRIALS) 
54 FORMAT (27HOLIMIT ON TRIALS EXCEEDED. ) 
55 FORMAT ('OBEST VERTEX IS NO.',I3,' OLD MIN WAS ' 
1E15.7, ' NEW MIN IS ',E15.7) 
56 FORMAT ('OMIN OBJECTIVE FUNCTION IS ',E15.7) 
END 
SUBROUTINE FBV (K,FUN,M) 
DIMENSION FUN(50) 
Meee 


> 


DO) 1 eZ 
IF (FUN(M)eLE FUNG)» GOuLomr 
M= TI 

1 CONTINUE 


RETURN 

END 

SUBROUTINE BOUT (NT,NPT,NFE,NCE,NV,NVT,V,K,FN,C, IK) 
DIMENSION V(50,50), FN(50), C(25) 

WRITE (6,4) NT,NPT,NFE,NCE 


DO 1 I=1,k 
WRITE (6,5) FN(L),(V(J,1),J=1,NV) 
IF (NVT.LE.NV) GO TO l 
NVP = NV+1 
WRITE (6,6) (V(J,1I),J=NVP,NVT) 
1 CONTINUE 


LF (RANE Cp Gem ene 


WRITE (6,7) (C(L),1I=1,NV) 
RETURN 


eZ 


2 IF (IK.GE.0) GO TO 3 
WRITE (6,8) (C(1L),I=1,NV) 


RETURN 

GREGG aL (Gye = 1 NV) 
RETURN 

4, FORMAT ('ONO. TOTAL TRIALS = ',15,4X,'NO. FEASIBLE 
ligase eee 

215,4X,'NO. FUNCTION EVALUATIONS = ',15,4X,'NO. 
3CONSTRAINT EVALUATIONS 

4= ',25/'0 FUNCTION VALUE’ ,6X, ' INDEPENDENT 

5 VARIABLES / DEPENDENT 


60R IMPLICIT CONSTRAINTS ' ) 

HORM AGU HS 57 52X,7E14.7/(21X,7E14.7)) 

FORMAT (21X,7E14.7) 
FORMAT (LOHOCENTROID 11X,7E14.7/(21X,7E14.7)) 
MOR an@en BEG VERTEX .7X,7E14.7/(21X,7E14.7)) 
HOA V@M@CENIROlD LESS VX’ ,12,2X%,)7£14.7/ (21X,7E£14.7)) 
END ; 

FUNCTION FE(X) 

DIMENSION X(3): 

COMMON TDIFF 

CALL PLANT(X) 

FE=TDIFF 

RETURN 

END 

FUNCTION KE(X) 

DIMENSION X(3) 

en) 

RETURN 

END 
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